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. ABSTRACT 

' We study non-linear axisymmetric pulsations of rotating relativistic stars using a gen- 

eral relativistic hydrodynamics code under the assumption of a conformally flat three- 
metric. We compare the results of our simulations, in which the spacetime dynamics 
' is coupled to the evolution of the fluid, to previous results performed in the Cowl- 

. ing approximation in which the spacetime dynamics was neglected. We show that 

the conformal flatness condition has only a small effect on the dynamics of pulsating 
relativistic stars and the obtained pulsation frequencies are very close to those ex- 
pected in full general relativity. The pulsations are studied along various sequences of 
both uniformly and differentially rotating relativistic polytropes with index iV = 1. 
^Y^' , For small pulsation amplitudes we identify several modes, including the lowest-order 

^ = 0, 2, and 4 axisymmetric modes, as well as several axisymmetric inertial modes. 
Differential rotation significantly shifts mode frequencies to smaller values, increasing 
lO i the likelihood of detection by current gravitational wave interferometric detectors. We 

' observe an extended avoided crossing between the I ~ and I = 4 first overtones 

I (previously known to exist from perturbative studies), which is important for cor- 

Q^. rectly identifying mode frequencies in case of detection. For uniformly rotating stars 

I ' near the mass-shedding limit, we confirm the existence of the mass-shedding- induced 

O I damping of pulsations and argue that it is still relevant for secularly unstable modes, 

even though the effect is not as strong as was previously found in the Cowling ap- 
^ ' proximation. We also investigate non-linear harmonics of the linear modes and notice 

that rotation changes the pulsation frequencies in a way that would allow for vari- 
ous parametric instabilities between two or three modes to take place. Although this 
scenario has been explored before for slowly-rotating collapse, it could become very 
interesting in the case of rapidly rotating collapse, where the quasi-radial mode could 
' be in resonance with inertial modes. We assess the detectability of each obtained mode 

by current gravitational wave detectors and outline how the empirical relations that 
have been constructed for gravitational wave asteroseismology could be extended to 
include the effects of rotation. 

Key words: Hydrodynamics - relativity - methods: numerical - stars: neutron - 
stars; oscillations - stars: rotation 
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1 INTRODUCTION 

Axisymmetric pulsations of rotating neutron stars could be 
excited in a number of astrophysical scenarios, such as ro- 
tating core collapse, accretion-induced collapse, core quakes 
due to a large phase-transition in the equation of state 
(EOS), or hypermassive neutron star formation in a binary 
neutron star merger (see Stergioulas, Apostolatos & Font 
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2004; Kokkotas & Stergioulas 2005, for extensive discus- 
sions) . These pulsations are a potential source of detectable 
high-frequency gravitational waves. While for nonrotating 
stars the frequencies of normal modes can be computed with 
perturbative methods and a theory of gravitational wave as- 
teroseismology has already been formulated (Andersson & 
Kokkotas 1998; Kokkotas, Apostolatos & Andersson 2001; 
Benhar, Ferrari & Gualtieri 2004), there exist no accurate 
frequency determinations for rapidly rotating stars to date, 
nor has the theory of gravitational wave asteroseismology 
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been extended to include the effects of rotation on the oscil- 
lation frequencies. Most existing computations of oscillation 
modes in rapidly rotating relativistic stars use the Cowl- 
ing approximation (Yoshida ot al. 2002; Stergioulas, Apos- 
tolatos & Font 2004; Yoshida, Yoshida & Eriguchi 2005), 
with the only exception being the computation of the two 
lowest-order quasi-radial modes in Font et al. (2002) (see 
also Shibata 2003; Shibata & Sekiguchi 2003). 

In this paper we present an accurate determination of 
several axisymmetric pulsation modes of rotating staxs in 
general relativity. The accurate knowledge of the frequen- 
cies of different modes excited in the astrophysical events 
mentioned above is necessary both for narrow-banding of 
the detectors, as well as for solving the inverse problem, i.e. 
identifying the EOS of high-density matter. 

The traditional approach for computing mode frequen- 
cies uses perturbation theory for either solving a time- 
independent eigenvalue problem or for obtaining the time 
evolution of the linearised equations governing the dynam- 
ics of matter and spacetime (see Kokkotas & Schmidt 1999; 
Kokkotas & Ruoff 2003, for comprehensive reviews of these 
approaches). The advantage of the pcrturbativc approach is 
that the equations can be expanded in terms of spherical 
harmonics. However, for rapidly rotating relativistic stars 
this approach has only worked in the Cowling approxima- 
tion so far (Yoshida & Eriguchi 1999, 2001; Yoshida et al. 
2002; Yoshida, Yoshida & Eriguchi 2005), except for zero- 
frequency /-modes, which were computed in full general rel- 
ativity (Stergioulas & Friedman 1998; Morsink, Stergioulas 
& Blattnig 1999). The main problem for applying the pcr- 
turbativc approach in full general relativity is the absence of 
analytic boundary conditions at infinity, which would allow 
to apply the outgoing-wave boundary conditions defining 
the quasi- normal modes (see Stergioulas 2003, for a review). 
Only if one assumes the slow-rotation approximation the 
problem is still tractable (see, e.g., Hartle & Friedman 1975; 
Kojima 1997; Datta et al. 1998; Ruoff & Kokkotas 2002; 
Pons et al. 2005). 

In recent years, the time evolution of the non-linear 
equations governing the dynamics of matter and space- 
time hcis been introduced as a promising new approach for 
computing mode frequencies (Font, Stergioulas & Kokko- 
tas 2000; Font et al. 2001, 2002; Stergioulas & Font 2001; 
Stergioulas, Apostolatos & Font 2004). For small ampli- 
tudes, the obtained frequencies are in excellent agreement 
with those expected by linear perturbation theory, while 
two-dimensional eigenfunctions can be obtained through a 
Fourier transform technique (see Stergioulas, Apostolatos & 
Font 2004, SAF hereafter). The advantages of this method 
are that one does not need precise outer boundary condi- 
tions and that one can also study non-linear pulsations. 

This study extends the results presented in SAF (which 
wore obtained in the Cowling approximation) by incorpo- 
rating the spacetime dynamics in the evolutions. This is 
done by using the Isenberg- Wilson-Mathews approximation 
of general relativity (also known as the conformal flatness 
condition; hereafter CFC) where the 3-1-1 Einstein equa- 
tions reduce to a non-linear sot of five coupled elliptic equa- 
tions for the lapse function, the shift vector, and the confor- 
mal factor (Isenberg 1978; Wilson, Mathews & Marronetti 
1996). The approximation essentially ignores gravitational 
radiation and is thus appropriate for equilibrium, quasi- 



equilibrium, but also for highly dynamical situations (see, 
e.g., Cook, Shapiro & Teukolsky 1996; Dimmelmeier, Font & 
Miiller 2002a; Oechshn, Rosswog & Thielemann 2002; Faber, 
Grandclement & Rasio 2004; Saijo 2004). 

For small pulsation amplitudes we identify several 
modes, including the lowest-order I = 0, 2, and 4 axisym- 
metric modes as well as several axisymmetric inertial modes. 
The pulsations are studied along the same sequences of uni- 
formly and differentially rotating relativistic polytropes with 
index A'^ = 1 as in SAF. Differential rotation significantly 
shifts mode frequencies to smaller values, increasing the like- 
lihood of detection by current gravitational wave interfer- 
ometric detectors. An important feature of the frequency 
spectrum, induced by rotation, is the existence of avoided 
crossings between different mode sequences (see Clement 
1986; Yoshida & Eriguchi 2001). We observe an extended 
avoided crossing between the / = and I = 4 first overtones. 
This is important for correctly identifying mode frequencies 
in case of detection. 

Our non-linear approach allows us to identify non-linear 
harmonics in addition to the well known linear modes. These 
harmonics arise due to couplings between various modes 
or due to non-linear self-couplings (see also Sperhake, Pa- 
padopoulos & Andersson 2001; Sperhake 2002); for similar 
results obtained for the non-linear oscillations of a torus or- 
biting a black hole, see Zanotti et al. (2005). It has been 
suggested (Clark 1979; Eardley 1983) that nonradial oscilla- 
tions after core bounce could be enhanced through a para- 
metric instability with the quasi-radial mode (see also Pas- 
samonti et al. 2005, for recent related work). In nonrotating 
or slowly rotating collapse, such a parametric instability can 
only take place under special conditions that would allow 
the two modes to be in resonance. In our work we find that 
rotational shifting of the frequency of different modes broad- 
ens the range of parameters for which interesting resonances 
could take place. In particular we notice that the quasi-radial 
mode will be in resonance with some inertial mode(s) for all 
rotation rates above a critical value. It is thus interesting to 
further study the possible energy transfer between different 
modes excited after, e.g., a core collapse or an accretion- 
induced collapse event, either on secular time-scales or as a 
possible parametric instability. 

The paper is organized as follows. In Section 2 we intro- 
duce the mathematical and numerical framework and specify 
the initial fluid perturbations, while in Section 3 we present 
the sequences of equilibrium initial models. In Section 4 we 
discuss the effects of linear pulsations, focusing on the role of 
rotation and avoided crossings of modes, and present a de- 
tailed comparison to results in the Cowling approximation. 
Section 5 is devoted to the technique of mode recycling, 
and in Section 6 we examine non-linear effects of the pul- 
sations like mode coupling or mass-shedding-induced damp- 
ing. Gravitational wave emission and asteroseismology are 
discussed in Section 7. A summary of our results in Section 8 
concludes this work. 

Unless otherwise noted, we choose dimensionless units 
for all physical quantities by setting the speed of light, the 
gravitational constant, and the solar mass to one, c = G = 
Mq = 1. Latin indices run from 1 to 3, Greek indices from 
1 to 4. 
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2 MATHEMATICAL AND NUMERICAL 
FRAMEWORK 

We study axisymmetric pulsations of rapidly rotating rela- 
tivistic stars by first constructing several sequences of uni- 
formly and differentially rotating equilibrium models. In 
SAF the equilibrium models were constructed using the nu- 
merical code rns (Stergioulas & Friedman 1995). In the 
present work we build the stellar equilibrium models us- 
ing the self-consistent field method described in Komatsu, 
Eriguclii & Hachisu (1989a, b) (KEH hereafter), which solves 
the general relativistic hydrostatic equations for rotating 
matter distributions whose pressure obeys an EOS given 
by a polytropic relation (see Eq. (9) below) . Comparisons of 
the accuracy of both approaches in the case of uniform rota- 
tion can bo found in Nozawa et al. (1998) and in Stergioulas 
(2003). Specific details of the equilibrium models are de- 
ferred to Section 3 below, where a quantitative comparison 
of the equilibrium properties of a particular highly differ- 
entially rotating model built either using rns or the KEH 
solver is made. 

These equilibrium models axe taken as initial data for 
the evolution code after a suitable perturbation has been 
added in order to excite specific modes of oscillation (sec 
Section 2.4). The time dependent numerical simulations are 
performed with the code CoCoNuT, developed by Dim- 
melmeier. Font & Miiller (2002a,b) with a metric solver 
based on spectral methods as described in Dimmelmeier 
et al. (2005). The code uses the general relativistic field 
equations for a curved spacetime in the 3 + 1-split under 
the assumption of conformal flatness for the three-metric. 
The hydrodynamics equations are consistently formulated in 
conservation form, and are solved by high-resolution shock- 
capturing schemes. 

In the code used by SAF, the spacetime dynamics was 
neglected, and the attention was focused on the oscillations 
of the stars on a fixed background metric given by the so- 
lution of the Einstein equations on the initial time slice. 
Keeping the spacetime fixed to the initial equilibrium state 
during the evolution corresponds to the Cowling approxi- 
mation in perturbation theory. In the simulations presented 
in this work, however, such a simplification is not made and 
the spacetime fields are also allowed to evolve in time. 

In the following, we present the mathematical formula^ 
tion of the metric and hydrodynamics equations, and then 
summarise the numerical methods used for solving them. 

2.1 Metric equations 

We adopt the ADM 3+1 formalism by Aruowitt, Deser & 
Misner (1962) to foliate a spacetime endowed with a metric 
g^v into a set of non-intersecting spacelike hypersurfaces. 
The line element then reads 

ds^ = gij,v dx'^dx^ = —a^ dt^ -\-^ij{dx'' + (f dt){dx^ + 13^ dt) , (1) 

where a is the lapse function, /3' is the spacelike shift three- 
vector, and is the spatial three-metric. 

In the 3 + 1 formalism, the Einstein equations are split 
into evolution equations for the three-metric 7ij and the ex- 
trinsic curvature Kij, and constraint equations (the Hamil- 
tonian and momentum constraints) which must be fulfilled 
at every spacelike hypersurface: 



dt'^ij = -2aKij + Vi/3j + Vj/3i, 

dtKij = -ViVja + a{Rij-2KikK^) + P^VkKij + KiuVjf}'' 
+KjkVip'' - 87ra(Sij --^{S^- Ph)) , (2) 
= R - Kij K'^ -167rfm, 

where the maximal slicing condition, Kl = 0, is imposed. In 
these equations Vi is the covariant derivative with respect to 
the three- metric 7;^, Rij is the corresponding Ricci tensor, 
and R is the scalar curvature. 

The matter fields of the general relativistic fluid ap- 
pearing in the above equations are the spatial compo- 
nents Sij of the (perfect fluid) stress-energy tensor T'"' = 
phu^u" + Pg'^" , the three momenta = phW'^v\ and the 
total energy pn = a'^T°°. The fluid is specifled by the rest- 
mass density p, the four- velocity , and the pressure P, 
with the speciflc enthalpy defined as /i = 1 + e + P/p, 
where e is the specific internal energy. The three-velocity 
of the fluid as measured by an Eulerian observer is given 
by = u''/{avP) + /3'/a, and the Lorentz factor W = au° 
satisfles the relation W = l/\/l — ViV'. 

The equations of the original ADM formulation, when 
implemented numerically with standard flnite-difference 
methods, suffer from several numerical instabilities. For 
many years there have been numerous attempts to refor- 
mulate these equations into forms better suited for numer- 
ical work (sec, e.g., Alcubierre et al. 2004, and references 
therein) . One recent approach is based on a constrained evo- 
lution scheme (Bonazzola et al. 2004), which exploits the 
observation that the more constraints are used in the for- 
mulation of the equations the more numerically stable the 
evolution appears to be. We refer the interested reader to 
the discussion in Dimmelmeier et al. (2005) and references 
therein. 

Based on the ideas of Isenberg (1978) and Wilson, 
Mathews & Marronetti (1996), and as it was done in the 
work of Dimmelmeier, Font & Miiller (2002a,b), we follow a 
similar strategy. We approximate the general metric p^j^ by 
replacing the spatial three-metric 7ij with the conformally 
flat three-metric 

7ij = 4>%j, (3) 

where 7i-j is the flat metric and </> is the conformal factor. 
Therefore, at all times during a numerical simulation we 
assume that all off-diagonal components of the three-metric 
are zero, and the diagonal elements have the common factor 

In this CFC approximation the expression for the ex- 
trinsic curvature becomes time-independent and reads 

Kij = ^ {Vipj + Vjpi - \iijVkfi''^ . (4) 

With this the ADM equations (2) reduce to a set of flve cou- 
pled elliptic non-linear equations for the metric components. 
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A<p = -2-k4>' phW' -P + 



r2 r, , K^,K'' 



167r 



U = [D,Si,T\, 



A{a(t>) = 27ra/ p/i(3iy^ - 2)+5P+ 



A/3' = 167ra0*5* + 2(j)^°K''Vj ( ^ ] - iv'V;t/3^ 



167r 

' ^4fiJ 3 



(5) 



where V, and A are the flat space Nabla and Laplace opera- 
tors, respectively. They do not contain explicit time deriva- 
tives, and thus the metric is calculated by a fully constrained 
approach, at the cost of neglecting some evolutionary de- 
grees of freedom in the spacetime metric (e.g., dynamical 
gravitational wave degrees of freedom). On each time slice 
the metric is hence solely determined by the instantaneous 
hydrodynamic state, i.e. the distribution of matter in space. 

The accuracy of the CFC approximation has been 
tested in various works, both in the context of stellar core 
collapse and for equilibrium models of neutron stars (Cook, 
Shapiro & Teukolsky 1996; Dimmelmeier, Font & Miiller 
2002a; Shibata & Sckiguclii 2004; Dimmelmeier ct al. 2005; 
Saijo 2005; Cerda-Duran et al. 2005). The spacetime of 
rapidly (uniformly or differentially) rotating neutron star 
models is still very well approximated by the CFC met- 
ric (3). The accuracy of the approximation is expected to 
degrade only in extreme cases, such as a rapidly rotating 
black hole. 

Recently, Cerda-Duran et al. (2005) have extended the 
CFC system of equations by incorporating additional de- 
grees of freedom in the approximation, which render the 
spacetime metric exact up to the second post-Newtonian or- 
der (CFC-f- approach). Results for uniformly rotating pulsat- 
ing neutron stars show only minute differences with respect 
to the CFC approximation for the computed frequencies of 
the quasi-radial fundamental F-mode and its first two over- 
tones. Moreover, a direct comparison of the CFC approach 
with fully general relativistic simulations of the quasi-radial 
modes of a rotating star, obtained in Font et al. (2002) have 
also yielded excellent agreement in the oscillation frequen- 



(8) 



o,r'"' 



)d\na 



respectively. Hero {)* = — P' /(-^, and g = 0^/7, with 
g = deX,{g^v) and 7 = det(7ij). In addition, F^^ are the 
Christoffel symbols associated with g'^,,. 

The system of hydrodynamics equations (7) is closed by 
an EOS, which relates the pressure to some thermodynam- 
ically independent quantities, e.g., P = P{p,e). For rotat- 
ing neutron star models below the mass-shedding limit (see 
Friedman, Ipser & Parker 1986, for the precise definition of 
the mass-shedding limit in the case of rapidly rotating rel- 
ativistic stars) we assume that the star remains isentropic. 
We can thus demand that the pressure obeys an EOS given 
by the polytropic relation 



P = Kp\ 



(9) 



where K is the polytropic constant and 7 = 1 -|- 1 /N is the 
adiabatic index. Note that the evolution equation for the 
generalised energy r can bo discarded if an EOS of the form 
P — P{p) as in Eq. (9) is used. In this particular case the 
internal specific energy can be obtained from the ideal fluid 
EOS as e = P/[p(7- 1)]. 

Near the mass-shedding limit, even small amplitude os- 
cillations can result in significant shedding of matter from 
the stellar surface in the form of shocks (see Section 6.3). 
In this case the polytropic relation (9) does not hold any- 
more. Therefore, as in the work by SAF, we then employ 
the adiabatic ideal fluid EOS instead. 



: pe(7 - 1). 



(10) 



We point out that at the initial time the isentropic equilib- 
rium models constructed with the polytropic EOS (9) are 
consistent with the ideal fluid EOS (10). 



2.2 General relativistic hydrodynamics 

The hydrodynamic evolution of a relativistic perfect fluid 
is determined by a system of local conservation equations, 
which read 



V^J" = 0, 



V^T"" = 0, 



(6) 



where J** = pu'^ is the rest-mass current, and V^j denotes 
the covariant derivative with respect to the four- metric g^^. 
Following Banyuls et al. (1997) we introduce a set of con- 
served variables in terms of the primitive (physical) variables 
(P, Vi,e): 



D = pW, Si = phW^Vi 



phW^ 



D. 



Using the above variables, the local conservation laws (6) 
can be written as a flrst-order, flux-conservative hyperbolic 
system of equations. 



dt 



(7) 



with the state vector, flux vector, and source vector 



2.3 Numerical methods for solving the metric and 
hydrodynamics equations 

The hydrodynamics solver performs the numerical time in- 
tegration of the system of conservation equations (7) using a 
high-resolution shock-capturing (HRSC) scheme on a finite 
difference grid (for a review of such methods in numerical 
general relativity, see Font 2003). This method ensures nu- 
merical conservation of physically conserved quantities and 
a correct treatment of discontinuities such as shocks (which 
may be present in hydrodynamic quantities). In (upwind) 
HRSC methods a Riemann problem has to be solved at each 
cell interface, which requires the reconstruction of the prim- 
itive variables {p,v',e) at these interfaces. We use the PPM 
method for the reconstruction, which yields third order accu- 
racy in space. The solution of the Riemann problems then 
provides the numerical fluxes at cell interfaces. To obtain 
this solution, the characteristic structure of the hydrody- 
namics equations is explicitly needed (Banyuls et al. 1997). 
In our code the numerical fluxes are computed by means of 
Marquina's approximate flux formula (Donat et al. 1998). 
The time update of the conserved vector U is done using 
the method of lines in combination with a Runge-Kutta 



© 2005 RAS, MNRAS 000, 1-24 



Pulsations of rotating relativistic stars 5 



scheme with second order accuracy in time. Once the state 
vector is updated in time, the primitive variables are recov- 
ered through an iterative Newton-Raphson method. 

Although the most common approach to numerically 
solve the Einstein equations is by means of finite difi^erences, 
such methods are not particularly well suited if the metric 
equations are formulated as non-linear coupled elliptic equa- 
tions like in the CFC approach. For multidimensional sim- 
ulations, the necessary grid resolutions typically require to 
solve computationally very expensive numerical problems. If 
iterative solvers in combination with spherical polar coordi- 
nates are used, an additional obstacle may manifest itself in 
slow or failed convergence due to problems at the coordinate 
origin or axis (see discussion in Dimmelmeier et al. 2005). 
A possible remedy to these shortcomings can be the use of 
non-linear multigrid solvers. 

As an alternative strategy to reduce the complexity 
associated with solving elliptic equations by reducing the 
number of grid points required for a given numerical accu- 
racy, we utilise am iterative non-linear solver based on spec- 
tral methods in our code. This metric solver, which is de- 
scribed in detail in Dimmelmeier ct al. (2005), uses routines 
from the publicly available object-oriented Lorene library 
(www.lorene.obspm.fr), which supplies routines that im- 
plement spectral methods in spherical polar coordinates. In 
contrast to the hydrodynamic quantities, the metric com- 
ponents are always smooth, and thus spectral methods are 
ideally suited for numerically representing the spacetime 
metric. The practical implementation of the combination of 
HRSC methods for the hydrodynamics and spectral meth- 
ods for the metric equations (the Mariage des Maillages or 
'grid wedding' approach) in a multidimensional numerical 
code has been presented in Dimmelmeier et al. (2005). 

The CoCoNuT code utilises Eulerian spherical polar 
coordinates {r,6,ip}, and thus axially or spherically sym- 
metric configurations can be easily simulated. For the rotat- 
ing neutron star models discussed in this work, we choose 
an axisymmetric grid setup (no dependence of quantities on 
the coordinate tp), and assume symmetry with respect to 
the equatorial plane. The finite difi'erenco grid consists of 
160 radial and 60 angular grid points, which are equidis- 
tantly spaced. A small part of the grid covers an artificial 
low-density atmosphere extending beyond the stellar sur- 
face, whose rest-mass density is 10^^^ of the initial central 
rest- mass density of the star. The spectral grid of the metric 
solver is split into 3 radial domains with 33 radial and 17 an- 
gular collocation points each. The innermost radial domain 
(or nucleus) stretches from the coordinate origin to half the 
stellar equatorial radius, followed by the second radial do- 
main which extends out to the outer boundary of the finite 
difference grid. The third domain uses a compactified radial 
coordinate and reaches out to spatial infinity. The metric 
equations (5) can therefore be numerically integrated out to 
spacelike infinity, and all noncompact support source terms 
in these equations can bo consistently handled in a non- 
approximative way. For further details about the grid setup 
and particularly the relevance of a compactified radial spec- 
tral grid for an accurate evolution of rotating equilibrium 
models, we refer to Dimmelmeier et al. (2005). 

Even when using spectral methods the calculation of 
the spacetime metric is computationally expensive. Hence, 
in our simulations the metric is updated only once every 



50 hydrodynamic time steps during evolution (which corre- 
sponds to a time interval of 10~^ ms) and extrapolated in 
between. The suitability of this procedure is tested and dis- 
cussed in detail in Dimmelmeier, Font & Miillcr (2002a) . Wo 
also note that convergence tests with different grid resolu- 
tion have been performed to ascertain that the regular grid 
resolution specified above is appropriate for our simulations. 

2.4 Fluid perturbations for exciting the pulsations 

To excite specific eigenmodes of oscillation in the stellar 
models we perturb selected equilibrium variables before 
starting the evolution. In the absence of the true eigenfunc- 
tion of a given mode, each perturbation is selected so as to 
mimic the angular dependence of the eigenfunction of the 
corresponding mode of a slowly rotating Newtonian star. 
Usually, this ensures that the chosen mode will dominate 
the time evolution at least for the slower rotating models. 
However, since the perturbation is not exact, additional pul- 
sation modes will be excited, especially for rapidly rotating 
models (see Section 5). 

For the / = modes we use as trial eigenfunction a 
perturbation of the radial component of the covariant three- 
velocity in the form 

«'--««-(-^). (11) 

where rg(^) is the coordinate radius of the surface of the 
star. The constant a is the amplitude of the perturbation, 
for which we choose the value —0.005 (in units of c). 

The I = 2 modes are excited by perturbing the 0- 
component of the covariant three-velocity as follows: 

ve = a sin I tt — ) sinS cos^, (12) 

where we set a to 0.01 (see Font et al. 2001, for more details 
about perturbations of this form). 

Correspondingly, the I — 4 perturbations are excited 
using the following perturbation: 

V0 = asm (^T^:^rj0yj ^'■^^ cos6'(3- 7cos^6') . (13) 

The time scries of the evolved perturbations are Fourier 
analysed after an evolution time of 20 ms, which corresponds 
to a nominal frequency resolution of 0.05 kHz. Due to the 
finite evolution time and the numerical viscosity which is 
damping the oscillations, the frequency peaks are not 5- 
functions, but correspond to a nearly symmetric bell shape, 
spreading over several frequency bins. As in SAF, we use a 
second-order accurate numerical derivative formula in order 
to find the location of the maximum of a frequency peak, 
when the peaks are well separated. In practice, this method 
gives agreement with expected frequencies which is signifi- 
cantly better than the nominal frequency resolution due to 
the finite evolution time. 

The peaJis in the Fourier spectra are identified with spe- 
cific pulsation modes, starting from the nonrotating member 
of the sequence, where the pulsation frequencies are known 
from perturbation theory. As the rotation rate increases, it 
becomes gradually more difficult to identify specific modes 
in the Fourier spectrum. For this reason, we also extract the 
two-dimensional eigenfunction for each peak in the Fourier 
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spectrum and use it as an additional criterion to identify 
specific modes (see Section 5). 

We also emphasize that in contrast to previous work 
(as, e.g., in Font ct al. 2001, and SAF), which assumed the 
Cowling approximation of a fixed spacetime metric, we do 
not use a perturbation of the rest-mass density for the I = 
mode, 

6p = apc sin (^7v^^^ , (14) 

where pc is the central rest-mass density of the star. In the 
Cowling approximation, the oscillation amplitude ttovoi of 
the density at a particular location in the star during the 
evolution not only scales linearly with the initial perturba- 
tion amplitude a (as expected for small o), but its value is 
also typically close to a. In contrast to this, if the spacetime 
is allowed to evolve, i.e. it is coupled to the evolution of the 
fluid, we find that the oscillation amplitude aovoi during the 
evolution is significantly larger than a, usually by a factor 
~ 5-8. Thus even a small initial perturbation with an ampli- 
tude of a few per cent can excite an oscillation with a large 
amplitude aevoi, which can possibly violate the assumption 
of linearity. 

We attribute this effect to the fact that the perturbar 
tion enters the hydrodynamic evolution equations (7) not 
only through the state vector U, but also via the metric 
in the form which is proportional to 6''. If U and (via 
the metric equations 5) also </)® are affected by the initial 
perturbation with amplitude a, the combination i/yC/ ex- 
hibits a much larger oscillation amplitude than a, namely 
flcvoi ~ (1 + a)^ • (1 + a) - 1 = (1 + af - 1, which is 
roughly 7a for a <C 1. This is also reflected by an increase 
of the total mass of the star by (1 -|- a)^, when applying this 
artificial perturbation by adding rest-mass density to the 
equilibrium model, as the rest mass integral also contains a 
term ^/7p. As the spacetime metric adapts to a new quasi- 
equilibrium state with this higher rest mass, the oscillation 
frequencies are systematically shifted towards higher values 
than those obtained if the approximately 'mass neutral' per- 
turbation (11) is applied. Note that this effect of frequency 
shift is negligible in the Cowling approximation when Sp is 
small, as there y'7 is constant in time, and thus the total 
rest mass of the star is affected by the density perturbation 
only through p, i.e. like {1 + a). 



3 EQUILIBRIUM MODELS 

Since our focus is on the effects of rotation on pulsation 
modes, we do not survey here a broad range of high-density 
EOSs, but choose instead a single polytropic EOS (9) with 
A'^ = 1 and K = 100. For the sake of comparison we use 
the same equilibrium model sequences as SAF to which we 
refer the reader for more details. Here we only give a brief 
overview of the basic properties of the various sequences. 

We restrict our attention to two different sequences of 
differentially rotating models (sequences A and B) and their 
uniformly rotating counterparts (sequences AU and BU). 
The equilibrium properties of all models are summarised 
in Table 1. In the nonrotating limit all sequences end in 
the same nonrotating model (thus, models AO, AUO, BO, 
and BUO all coincide). Notice that in Table 1 the numeri- 



cal values of different equilibrium properties are displayed 
as computed by SAF, using the numerical code rns (Ster- 
gioulas & Friedman 1995). In contrast, the equilibrium ini- 
tial data used in the actual time evolutions presented here 
are obtained by a similar, but different, initial data solver, 
based on the original KEH method (Komatsu, Eriguchi & 
Hachisu 1989a). An extensive comparison between the two 
initial data solvers shows very good agreement with differ- 
ences being much less than the 1 per cent level in all com- 
puted equilibrium quantities (the differences arise solely due 
to the truncated computational domain of the original KEH 
method). We also note that for a few of the most rapidly 
differentially rotating equilibrium models used by SAF, the 
elliptic solver for the spacetime evolution used here did not 
converge. These models are therefore omitted in Table 1. 

The differentially rotating sequence A and its corre- 
sponding uniformly rotating sequence AU are characterised 
by a fixed rest mass Mo — 1.506 M0. Along sequence A, 
the degree of differential rotation is held fixed ai A — 1, 
where A = A/ re with re being the equatorial coordinate 
radius of the star and A being the rotation parameter as 
defined in Komatsu, Eriguchi & Hachisu (1989a). The val- 
ues of Mo and A are chosen in order to represent a newly- 
born, differentially rotating neutron star. The angular ve- 
locity at the equator is roughly 1/3 to 1/2 of the central 
angular velocity, which is similar to the degree of differen- 
tial rotation obtained in typical core collapse simulations 
(see, e.g., Villain et al. 2004). The fastest rotating model in 
sequence A has a ratio of polar to equatorial coordinate axis 
of only rp/vc = 0.370, which translates into a rotation rate 
T/\W\ = 0.223. (Here T is the rotational kinetic energy and 
\W\ is the gravitational potential energy.) The central rest- 
mass density is nearly an order of magnitude smaller than 
the one of the corresponding nonrotating model, while the 
circumferential radius is nearly twice as large. The uniformly 
rotating sequence AU only reaches an axis ratio of 0.575 at 
a rotation rate r/|W^| = 0.095, half the central rest-mass 
density, and a 50 per cent larger radius than the correspond- 
ing nonrotating model. Model AU5 is at the mass-shedding 
limit . 

On the other hand, the differentially rotating se- 
quence B and its corresponding uniformly rotating se- 
quence BU are characterised by a fixed central rest-mass 
density pc = 1.28 x 10~^. In sequence B the differential ro- 
tation parameter is also set to A = 1. Its fastest rotating 
member has a gravitational mass M = 2.02 Mq and an axis 
ratio of 0.55. Since all models in the sequence are compact, 
the radius R is much smaller than along sequence A. The 
corresponding uniformly rotating sequence BU only reaches 
an axis ratio of 0.58 at the mass-shedding limit (model BU9) 
with an increase in radius R by 40 per cent. Thus, we see that 
when considering a sequence of fixed central rest-mass den- 
sity, the uniformly rotating models attain a larger equatorial 
radius than differentially rotating models, with the latter ex- 
panding out of the equatorial plane, becoming torus-like. 



4 LINEAR PULSATION MODES 

When the pulsation amplitude is small (so that, e.g., den- 
sity variations are at the level of Sp/p ~ 10~^), the dynamics 
is close to linear and one can identify the eigenfrequencies 
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Table 1. Properties of the four sequences of equilibrium models. A is a sequence of fixed rest mass Mq = 
1.506 Mq with A = 1, AU is the corresponding sequence of uniformly rotating models, B is a sequence of 
fixed central rest-mass density pc = 1.28 X 10~^ with A = 1, and BU is the corresponding sequence of 

uniformly rotating models. All models arc relativistic polytropes with = 1 and K = 100. Sc is the central 
energy density with s = p(l + e), M is the gravitational mass, -R is the circumferential stellar radius, and 
^c/e the central/equatorial angular velocity. All other quantities are defined in the main text. 



Model £c M R re Vp/re He f^e T/\W\ 

(xlO-3) (xlO-2) (xlO-2) 



AO 


1 


.444 


1.400 


9.59 


8.13 


1.000 


0.000 


0.000 


0.000 


Al 


1 


.300 


1.405 


10.01 


8.54 


0.930 


2.019 


0.759 


0.018 


A2 


1 


.187 


1.408 


10.40 


8.92 


0.875 


2.580 


0.977 


0.033 


A3 


1 


.074 


1.410 


10.84 


9.35 


0.820 


2.944 


1.125 


0.049 


A4 





.961 


1.413 


11.37 


9.87 


0.762 


3.192 


1.232 


0.066 


A5 





.848 


1.418 


12.01 


10.49 


0.703 


3.340 


1.303 


0.086 


A6 





.735 


1.422 


12.78 


11.25 


0.643 


3.383 


1.336 


0.107 


A7 





.622 


1.427 


13.75 


12.21 


0.579 


3.339 


1.337 


0.131 


A8 





.509 


1.433 


15.01 


13.45 


0.513 


3.197 


1.300 


0.158 


A9 


0, 


.396 


1.439 


16.70 


15.13 


0.444 


2.953 


1.223 


0.189 


AlO 





.283 


1.447 


19.03 


17.44 


0.370 


2.604 


1.101 


0.223 


AUO 


1.444 


1.400 


9.59 


8.13 


1.000 


0.000 


0.000 


0.000 


AUl 


1.300 


1.404 


10.19 


8.71 


0.919 


1.293 


1.293 


0.020 


AU2 


1.187 


1.407 


10.79 


9.30 


0.852 


1.656 


1.656 


0.037 


AU3 


1.074 


1.411 


11.56 


10.06 


0.780 


1.888 


1.888 


0.055 


AU4 


0.961 


1.415 


12.65 


11.14 


0.698 


2.029 


2.029 


0.076 


AU5 


0.863 


1.420 


14.94 


13.43 


0.575 


2.084 


2.084 


0.095 


BO 


1 


.444 


1.400 


9.59 


8.13 


1.000 


0.000 


0.000 


0.000 


Bl 


1 


.444 


1.437 


9.75 


8.24 


0.950 


1.801 


0.666 


0.013 


B2 


1, 


.444 


1.478 


9.92 


8.36 


0.900 


2.574 


0.944 


0.026 


B3 


1 


.444 


1.525 


10.11 


8.49 


0.850 


3.189 


1.160 


0.040 


B4 


1 


.444 


1.578 


10.31 


8.63 


0.800 


3.728 


1.342 


0.055 


B5 


1, 


.444 


1.640 


10.53 


8.77 


0.750 


4.227 


1.504 


0.071 


B6 


1 


.444 


1.713 


10.76 


8.91 


0.700 


4.707 


1.651 


0.087 


B7 


1 


.444 


1.798 


11.01 


9.05 


0.650 


5.185 


1.789 


0.105 


B8 


1 


.444 


1.899 


11.26 


9.17 


0.600 


5.683 


1.921 


0.124 


B9 


1 


.444 


2.020 


11.50 


9.26 


0.550 


6.232 


2.052 


0.144 


BUO 


1 


.444 


1.400 


9.59 


8.13 


1.000 


0.000 


0.000 


0.000 


BUI 


1 


.444 


1.432 


9.83 


8.33 


0.950 


1.075 


1.075 


0.012 


BU2 


1 


.444 


1.466 


10.11 


8.58 


0.900 


1.509 


1.509 


0.024 


BU3 


1 


.444 


1.503 


10.42 


8.82 


0.850 


1.829 


1.829 


0.037 


BU4 


1 


.444 


1.543 


10.78 


9.13 


0.800 


2.084 


2.084 


0.050 


BUS 


1 


.444 


1.585 


11.20 


9.50 


0.750 


2.290 


2.290 


0.062 


BU6 


1 


.444 


1.627 


11.69 


9.95 


0.700 


2.452 


2.452 


0.074 


BUT 


1, 


.444 


1.666 


12.30 


10.51 


0.650 


2.569 


2.569 


0.084 


BU8 


1 


.444 


1.692 


13.07 


11.26 


0.600 


2.633 


2.633 


0.091 


BU9 


1 


.444 


1.695 


13.44 


11.63 


0.580 


2.642 


2.642 


0.092 



and eigenfunctions of linear quasi-normal modes from a time 
evolution. In order to compute the real part of the eigenfre- 
quency of a pulsation mode, we Fourier-transform the time 
series of the evolution of a suitable physical variable (the 
density for the I = modes and ve for the I = 2 modes). 
Instead of examining the Fourier spectra at a few specific 
points inside the star, we integrate the amplitude of the 
Fourier transform along a coordinate line. For instance, in 
the case of / = modes we examine the integrated Fourier 
amplitude along 6 = Tr/2 (equatorial plane), while for the 



I = 2 modes the integrated Fourier amplitude along a line of 
6 = 7r/4 is used. We also verify that each identified discrete 
mode has the same frequency at any point inside the star in 
the coordinate frame. 

As discussed in SAF, the trial eigenfunction used for 
exciting the pulsations docs not correspond exactly to a 
particular mode. Therefore, additional modes apart from 
the main mode one wishes to study are excited, particu- 
larly for rapidly rotating models where rotational coupling 
effects are significant and higher-order coupling terms in the 
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Table 2. Frequencies of the fundamental quasi-radial {I = 0) 

mode, F, its first overtone, Hi , the fundamental quadrupole [l = 
2) mode, '^f, its first overtone, ^pi, and three inertial modes, i—2, 
ii , and 12 for the sequence of differentially rotating models A. All 
frequencies are given in kHz. 



Model 


F 


Hi 




2pi 


i-2 


ii 


«2 


AO 


1.458 


3.971 


1.586 


3.726 


0.000 


0.000 


0.000 


Al 


1.400 


3.816 


1.577 


3.580 


0.302 


0.460 


0.596 


A2 


1.358 


3.733 


1.567 


3.424 


0.399 


0.603 


0.779 


A3 


1.307 


3.664 


1.550 


3.237 


0.477 


0.711 


0.917 


A4 


1.248 


3.583 


1.535 


3.013 


0.543 


0.794 


1.022 


A5 


1.184 


3.494 


1.513 


2.780 


0.603 


0.863 


1.108 


A6 


1.105 


3.352 


1.482 


2.557 


0.646 


0.914 


1.163 


A7 


1.018 


3.360 


1.432 


2.315 


0.683 


0.946 


1.198 


A8 


0.915 


3.114 


1.360 


2.068 


0.711 


0.965 


1.204 


A9 


0.809 


2.985 


1.264 


1.827 


0.723 


0.958 


1.180 


AlO 


0.685 


2.830 


1.098 


1.610 


0.712 


0.911 


1.104 



Table 3. Same as Table 2, but for the sequence of uniformly 
rotating models AU. 



Model 


r 




-I 


"/'I 


/ 2 


'i 


h 


AUO 


1.458 


3.971 


1.586 


3.726 


0.000 


0.000 


0.000 


AUl 


1.398 


3.785 


1.562 


3.455 


0.280 


0.354 


0.468 


AU2 


1.345 


3.716 


1.554 


3.192 


0.384 


0.478 


0.611 


AU3 


1.283 


3.635 


1.537 


2.885 


0.468 


0.575 


0.736 


AU4 


1.196 


3.552 


1.516 


2.520 


0.545 


0.660 


0.810 


AU5 


1.107 


3.457 


1.459 


2.090 


0.639 


0.747 


0.858 



mode-eigenfunctions become comparable to the dominant 
term. We begin the identification of specific modes along 
a sequence of equilibrium models using the known pulsa- 
tion frequencies of the nonrotating star in the sequence and 
by comparing the Fourier transforms of evolved variables 
between subsequent models. Erroneous mode identification 
could happen close to the mass-shedding limit due to the ap- 
pearance of avoided crossings between different modes (sec 
Section 4.2). To avoid this we do not rely on Fourier trans- 
forms at only a few points inside the star, but in many cases 
reconstruct the whole two-dimensional eigcnfunction of each 
mode, using Fourier transforms at every point inside the 
star. At the cigenfrequency of a specific mode, the amplitude 
of the Fourier transform correlates with its eigcnfunction. 
A change in sign in the eigcnfunction corresponds to both 
the real and imaginary part of the Fourier transform going 
through zero. Comparing the eigenfunctions corresponding 
to different peaks in the Fourier transforms allows for an 
unambiguous identification of specific mode sequences. 

Tables 2 to 6 summarise our main results, showing the 
frequencies of the two lowest-order / — modes (F and 
Hi) and / = 2 modes (^f and ^pi) for all four sequences as 
well as the two lowest order I = 4 modes (*/ and ''pi) for 
sequence A. In addition, frequencies of several axisymmetric 
inertial modes are displayed. All of the above frequency data 
(except for the ''/-mode, which is omitted for reasons of 
clarity) are shown as a function of T/IVFI in Figs. 1 and 2. 
Next, we first discuss some general trends due to rotation 
and then present our results for each mode sequence in more 
detail. 



Table 4. Same as Table 2, but for the sequence of differentially 
rotating models B. 



i\ it X L(* 


1 ./ 


-"i 




in 


/ 2 


'i 


/ 2 


BO 


1.458 


3.971 


1.586 


3.726 


0.000 


0.000 


0.000 


Bl 


1.407 


3.927 


1.628 


3.713 


0.258 


0.399 


0.518 


B2 


1.373 


3.927 


1.670 


3.666 


0.370 


0.567 


0.739 


B3 


1.332 


3.964 


1.709 


3.584 


0.463 


0.702 


0.912 


B4 


1.287 


4.014 


1.747 


3.490 


0.544 


0.819 


1.061 


B5 


1.237 


4.072 


1.789 


3.390 


0.622 


0.921 


1.208 


B6 


1.178 


4.118 


1.819 


3.280 


0.681 


1.015 


1.318 


B7 


1.077 


4.179 


1.854 


3.187 


0.750 


1.090 


1.439 


B8 


1.080 


4.212 


1.900 


3.103 


0.827 


1.213 


1.562 


B9 


0.945 


4.469 


1.917 


3.028 


0.890 


1.294 


1.670 



Table 5. Same as Table 2, but for the sequence of uniformly 
rotating models BU. 



Model 


F 


Hi 


^/ 




i-2 


n 


«2 


BUO 


1.458 


3.971 


1.586 


3.726 


0.000 


0.000 


0.000 


BUI 


1.413 


3.915 


1.611 


3.634 


0.224 


0.288 


0.385 


BU2 


1.380 


3.907 


1.635 


3.516 


0.326 


0.413 


0.541 


BU3 


1.343 


3.921 


1.669 


3.345 


0.408 


0.511 


0.658 


BU4 


1.304 


3.950 


1.698 


3.200 


0.486 


0.601 


0.762 


BU5 


1.281 


3.964 


1.714 


3.018 


0.552 


0.675 


0.847 


BU6 


1.219 


4.010 


1.729 


2.859 


0.617 


0.743 


0.912 


BU7 


1.207 


4.018 


1.720 


2.677 


0.680 


0.809 


0.969 


BU8 


1.168 


4.030 


1.685 


2.512 


0.723 


0.868 


1.004 


BU9 


1.169 


4.029 


1.679 


2.483 


0.737 


0.878 


1.004 



4.1 General trends of rotational effects 

It is well known that the frequencies of the fundamental 
I = and / = 2 polar modes of oscillation depend mainly on 
the central density of a star, or, equivalently, on the com- 
pactness M/R (see, e.g., Hartle & Friedman 1975). The se- 
quences A and AU of fixed rest mass Mo — 1.506 Mq start 
with a nonrotating model with compactness M/R = 0.15 
and terminate at models with much smaller compactness 
[M/R = 0.076 for sequence A and M/R = 0.095 for s^ 
quence AU). Based on this significant decrease of the com- 
pactness along the fixed-rest-mass sequences, one expects a 

Table 6. Frequencies of the fundamental hexadecapole {I = 4) 
mode, and its first overtone, for the sequence of differen- 
tially rotating models A (left) and B (right). All frequencies are 
given in kHz. 



Model 






Model 






AO 


2.440 


4.896 


BO 


2.440 


4.896 


Al 


2.370 


4.727 


Bl 


2.453 


4.877 


A2 


2.300 


4.600 


B2 


2.468 


4.855 


A3 


2.223 


4.372 


B3 


2.486 


4.827 


A4 


2.130 


4.130 


B4 


2.500 


4.781 


A5 


2.028 


3.864 


B5 


2.504 


4.740 


A6 


1.910 


3.617 


B6 


2.491 


4.643 


A7 


1.780 


3.116 


B7 


2.499 


4.556 


A8 


1.630 


2.790 


B8 


2.501 


4.506 


A9 


1.480 


2.430 


B9 


2.493 


4.164 


AlO 


1.330 


2.028 
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Figure 1. Frequencies of various modes for sequences A (solid 
lines) and AU (dashed lines). Note the avoided crossing between 
the Hi and the *pi-mode. 



corresponding decrease in the frequencies of the fundamen- 
tal modes (and a similar tendency for the first overtones). 
In contrast, along sequences B and BU, where the central 
density is fixed, the compactness varies much less than for se- 
quences A and AU. In fact, for sequence B, the compactness 
even somewhat increases. One therefore expects a weaker 
dependence of the pulsation frequencies on rotation for the 
sequences of fixed central density. The above expectations 
have already been verified qualitatively in the Cowling ap- 
proximation by SAF. 

In slowly rotating stars, the frequencies of all inertial 
modes increase linearly with increasing r/| W| (see Friedman 
& Lockitch 2001; Lockitch, Friedman & Andersson 2003). 
At higher rotation rates, higher-order rotational terms can 
modify this behavior. For uniformly rotating stars, the ex- 
pectation is that the inertial mode frequency still increases 
up to the mass-shedding limit. As clearly visible for the 
rapidly rotating models of sequence A in Fig. 1, this gen- 
eral expectation is no longer valid for differentially rotating 
stars (for a more detailed discussion, see Section 4.4). 

Due to differential rotation, the outer layers of the star 
rotate slower and the equatorial radius is smaller compared 
to a uniformly rotating model of same r/|Vl^|. This leads to 
a smaller sound-crossing time and correspondingly higher 
fundamental mode frequencies for the differentially rotating 
models. This explains why the curves for the fundamental 
mode frequencies of the F and ^f-mode of sequence A in 
Fig. 1 have smaller slopes than those corresponding to se- 
quence AU. This behavior was already found for the same 
model sequence in the Cowling approximation in the work 
of SAF. 

In general, higher order or large I modes are affected 
more strongly by rotation than lower order or smaller I 
modes. At large rotation rates this can lead to avoided cross- 
ings between mode sequences, where modes can exchange 
the character of their eigenfunctions. These avoided cross- 



Figure 2. Same as Fig. 1, but for the sequences B and BU. 

ings are already known to exist from perturbative studies of 
axisymmetric modes in rotating Newtonian stars (Clement 
1986) and for relativistic quasi-radial modes (Yoshida & 
Eriguchi 2001). It is not trivial to decide how to label the 
mode sequences after an avoided crossing. The decisive cri- 
terion for labeling a pulsation mode is not the continuity of 
the eigenfrequency along a mode sequence. More important 
is the character of the oscillation, i.e. the eigenfunction. One 
must therefore examine the eigenfunctions of two pulsations 
before and after an avoided crossing. Then the character of 
the modes after the crossing can be determined according to 
which modes (of those before the crossing) they resemble. 
Thus, continuity of eigenfunctions is preferred over continu- 
ity of eigenfrequencies in labeling modes. 

4.2 Quasi-radial (Z = 0) modes and avoided 
crossings 

The computed frequencies for the fundamental quasi-radial 
/ = mode F and its first overtone Hi for the fixed rest mass 
sequences A and AU are displayed in Tables 2 and 3 and 
plotted in Fig. 1. Along the uniformly rotating sequence AU, 
there is a decrease in the frequency of the fundamental quasi- 
radial mode, since the central density of the star decreases 
with increasing rotation rate. For Newtonian nonrotating 
polytropic models the frequency of the fundamental quasi- 
radial mode is proportional to the square root of the average 
density. Even though we do not compute an average density 
for the rapidly rotating models of sequences A and AU, we 
notice that the decrease in the frequency of the _F-mode from 
1.458 kHz to 1.107 kHz along sequence AU follows closely 
the decrease in central energy density Sc- Along the differ- 
entially rotating sequence A, the frequency of the _F-mode 
is further decreasing, reaching a very low value of 685 Hz 
for the most rapidly rotating model. Remarkably, when one 
compares models of same r/|VK| along the two sequences A 
and AU, the frequency of the _F-mode is insensitive to the 
degree of differential rotation. 
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The first overtone Hi also decreases along the uni- 
formly rotating sequence AU, but by less than the near 
Y^-dependence of the F-mode. Along the differentially ro- 
tating sequence A, the iifi-modc shows a similar inscnsi- 
tivity to the degree of differential rotation as the F-mode, 
when comparing models of same T/|W^|, up to rotation rates 
of T/IM^I ~ 0.1. For larger rotation rates, however, an ex- 
tended avoided crossing with the *pi-mode takes place. Such 
avoided crossings exist, because the various mode sequences 
are affected to a different degree by rotation. In partic- 
ular, the frequencies of higher I modes tend to decrease 
faster with increasing T/|VF1 than the frequencies of lower 
I modes, leading to approaching mode-sequences. However, 
mode-sequences that contain similar terms in their eigen- 
function expansions are not allowed to cross, even in the 
linear approximation (Clement 1986; Yoshida & Eriguchi 
2001). Instead, two continuous sequences of pulsation fre- 
quencies avoid to cross, as shown in Fig. 1. 



The character of the eigenfunctions along these contin- 
uous sequences is exchanged at an avoided crossing. As a 
result, the lower- frequency part of the continuous sequence 
starting as the Hi-mode in the nonrotating limit becomes 
the *pi-modc at large T/|VF|. Correspondingly, the lower- 
frequency part of the continuous sequence which starts as 
the ^pi-modo iu the nonrotating limit, becomes the Hi- 
mode at large r/|VF|. At the avoided crossing, the char- 
acter of the eigenfunctions of both modes is a mixture of 
the eigenfunctions of the two modes before the avoided 
crossing. We have determined the correct labeling of the 
sequences after the avoided crossing by carefully compar- 
ing their two-dimensional eigenfunctions (see Section 5 for a 
discussion of how the eigenfunctions are obtained from our 
time-evolutions). The particular avoided crossing between 
the Hi and *pi-modcs was also observed along a sequence 
of rotating models in the relativistic Cowling approximation 
(Yoshida & Eriguchi 2001). 

Along the fixed central density sequences B and BU (see 
Tables 4 and 5 and Fig. 2), the frequency of the F-mode de- 
creases, because, even though the central density stays fixed, 
rotational effects still increase the radius of the star, so that 
the sound crossing time in the equatorial region increases. 
This seems to have a monotonous influence on the frequency 
of the F-mode, which is again insensitive to the degree of dif- 
ferential rotation. The frequency of the _ffi-ovcrtone, on the 
other hand, first decreases somewhat with rotation, but then 
increases again, eventually surpassing its value of 3.971 kHz 
in the nonrotating limit. Since the i?i-mode has a node in 
its eigenfunction in the nonrotating limit, it is more sensi- 
tive to rotational effects (in comparison to the fundamental 
-F-mode). These rotational effects have a different influence 
near the symmetry axis than near the equator. The chang- 
ing dependence of the i/i-mode frequency as the star be- 
comes more flattened reflects this sensitivity. In the case of 
sequence B, the avoided crossing between the _ffi-mode and 
the ''pi -mode happens at the largest rotation rates along this 
sequence. Consequently, the fastest rotating model B9 is in 
the avoided crossing region, where the mode eigenfunctions 
are mixed. 



4.3 Quadrupole (Z — 2) modes 

The computed frequencies for the fundamental quadrupole 
I = 2 mode and its first overtone ^pi for the fixed 
rest mass sequences AU and A are displayed in Tables 3 
and 2 and plotted in Fig. 1. Along the uniformly rotating 
sequence AU, there is only a small decrease in the frequency 
of the fundamental ^f-mode, from 1.586 kHz to 1.459 kHz 
at the mass-shedding limit. However, for the differentially 
rotating sequence A, the rate of decrease gets stronger for 
models with very high T/|VKj, and the frequency of the fun- 
damental ■^/-mode becomes as small as 1.098 kHz. 

The first overtone ^pi shows a much stronger decrease 
in frequency with increasing rotation rates. Starting from 
3.726 kHz at the nonrotating limit, along sequence AU its 
frequency becomes 2.090 kHz at mass-shedding, while it has 
a value of only 1.610 kHz for the fastest differentially rotat- 
ing model along sequence A. A striking difference compared 
to all other modes studied here is that the frequency of the 
^pi-mode is indeed sensitive to the degree of differential ro- 
tation, as is evident from Figs. 1 and 2. 

Along the fixed central density sequences B and BU 
the ^/-mode shows the opposite tendency compared to its 
behavior along the sequences A and AU, with its frequency 
increasing to 1.679 kHz and 1.917 kHz, respectively (see Ta- 
bles 4 and 5 and Fig. 2). In contrast, the frequency of the ^pi- 
mode is still decreasing (although not as drastically as for se- 
quences A and AU) with increasing rotation rates, reaching 
frequencies of 3.028 kHz at mass-shedding along sequence B 
and 2.483 kHz along sequence BU. The '^pi-modc appears to 
be even more sensitive to the degree of differential rotation 
along the fixed central density sequence B than along the 
fixed rest mass sequence A. 



4.4 Inertial modes 

In our simulations we observe a large number of inertial 
modes, which axe supported by the Coriolis force and be- 
come degenerate at zero frequency in nonrotating stars. 

Fi'oiii linear perturbation theory we expect that there exists 
an infinite number of inertial modes in a finite frequency 
range, which corresponds to to 2 f2 for uniformly rotat- 
ing Newtonian stars (see, e.g., Lockitch & Friedman 1999). 
In spite of this, we do not find evidence for the excitation 
of an arbitrary number of inertial modes, but only a few 
specific inertial modes are predominantly excited. Since we 
do not make use of the eigenfunction recycling method to 
excite inertial modes (see Section 5) , they are excited as by- 
products of the excitation of other modes. Hence, inertial 
modes can be excited either due to the non-exact nature 
of the initial trial eigenfunctions used to perturb the initial 
model (as described in Section 2.4) or due to non-linear cou- 
plings with linear polar modes (see discussion in Section 6.2). 
The fact that the amplitude of the observed inertial modes, 
with our choice of trial eigenfunctions, scales non-linearly 
with increasing amplitude, points to non-linear couplings as 
a possible origin, at least for most of the inertial modes we 
observe. 

Examining the eigenfunctions of the excited inertial 
modes we notice that the mode with the usually highest 
PSD in a Fourier transform has only one node in the eigen- 
function for ve, while other modes at larger or smaller fre- 
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Figure 3. Comparison of the frequencies for the ^/ and '^pi- 
mode (upper panel) and for the F and _ffi-mode (lower panel) 
in CFC (solid lines) and in the Cowling approximation (dashed 
lines) for the sequence of differentially rotating models A. In the 
lower panel we also show the frequencies of the '*pi-mode after 
the avoided crossing with the Hi-mode (green line). 



quency than this 'fundamental' inertial mode have a larger 
number of nodes (increasing as the absolute frequency dif- 
ference from the 'fundamental' inertial mode increases) . Out 
of the many excited inertial modes, we choose to display in 
Tables 2 to 5 and Figs. 1 and 2 the 'fundamental' inertial 
mode with only one node in ve, which we call ii-mode, and 
the two modes with two nodes in the eigenfunction of ve on 
both sides of the zi-mode in the frequency domain, which 
we denote the i_2-mode and i2-mode, respectively. 

Notice that the number of nodes is determined at slow 
rotation rates and can change at very high rotation rates. 
Note also that inertial modes come in two flavours (polar- 
led and axial-led) and are hybrid modes in the sense that 
they do not reduce to either a pure polar or a pure axial 
mode in the nonrotating limit (see Lockitch & Friedman 
1999). A proper classification scheme for inertial modes, 
based on the dominant components in their eigenfunctions, 
is already in use (Lockitch & Friedman 1999; Friedman 
& Lockitch 2001; Lockitch, Friedman & Andersson 2003). 
However, comparing our numerically calculated eigenfunc- 
tions with eigenfunctions obtained from linear perturbation 
theory with more simplifying assumptions (weak gravity or 
uniform density) is beyond the scope of the present paper 
and will be addressed in a separate publication. Thus, we 
use our own naming scheme for inertial modes simply for 
convenience, until proper identification is established in the 
future. 

Along sequence AU the frequencies of the three iner- 
tial modes increases monotonically with rotation. However, 
along sequence A the frequencies of the three inertial modes 
reach a maximum value for the rotation parameter r/|tV| 
being in the range 0.15 to 0.19, depending on the specific 
mode. This is due to the fact that, even though r/|H^| in- 
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Figure 4. Same as Fig. 3, but for the sequence of uniformly 
rotating models AU. 



creases monotonically along sequence A, both the angular 
velocity at the centre and the angular velocity at the sur- 
face reach maximum values along this sequence, as can be 
seen in Table 1. On the other hand, along sequences BU 
and B the frequencies of the three inertial modes again in- 
crease monotonically, as does the angular velocity at the 
centre and at the surface. 



4.5 Comparison to the Cowling approximation 

When comparing the dependence of the mode frequency on 
T/\W\ plotted in Figs. 3 to 6 for the F, H, V, and ^pi- 
modes obtained with the CFC approximation for the space- 
time evolution against previous results in the Cowling ap- 
proximation by SAF, several qualitative observations can 
be made. Apparently, for all sequences there is a large dis- 
crepancy in the values of the frequencies of the fundamen- 
tal quasi-radial _F-mode in the CFC simulations compared 
to the ones in the Cowling approximation at all rotation 
rates. This is also quantitatively evident in the relative dif- 
ferences of frequencies between the CFC and Cowling simu- 
lations of up to a factor ~ 2, listed for the nonrotating model 
(AO/AUO/BO/BUO) and the most rapidly rotating model of 
each sequence (AlO, AU5, B9, BU9) in Table 7. We notice 
that these findings are consistent with those of Yoshida & 
Kojima (1997) for nonrotating stars. In terms of absolute 
differences, for the _F-mode the Cowling approximation fails 
by far to predict the correct mode frequencies, as it yields 
values which are typically > 1 kHz too high. While in se- 
quences A, AU, and BU the frequency curve in the Cowling 
approximation at least correctly captures the decline with 
increasing rotation rate, it cannot reproduce the similar be- 
havior in sequence B. 

We emphasize that neither in sequence A nor B we find 
any evidence for a splitting of the F-mode, a phenomenon 
that was noticed for those sequences in the Cowling approxi- 
mation by SAF. This apparently confirms the possible expla- 
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Figure 5. Same as Fig. 3, but for the sequence of differentially 
rotating models B. 



Figure 6. Same as Fig. 3, but for the sequence of uniformly 
rotating models BU. 



nation for the _F-mode splitting as an artifact of the Cowling 
approximation offered in that work. In Figs. 3 to 6 and Ta- 
ble 7 we use the frequencies of the regular f-mode rather 
than those of the _/<ii-mode of SAF, as the eigenfunctions of 
the latter mode do not possess the typical characteristics of 
a genuine _F-mode (and, in addition, their amplitude of the 
frequency peak in FFTs is smaller than for the _F-mode). 

Figs. 3 to 6 demonstrate that for the frequencies of 
the fundamental quadrupole ^/-mode and the H\ and ^pi- 
modes, the discrepancies between the CFC simulations and 
the Cowling approximation are much less severe (see also 
Table 7). In the case of the _ffi-mode and the ^pi-mode, 
the curve in the Cowling approximation follows our CFC re- 
sults reasonably close, while we observe a crossing of the 
frequency curves for the ^/-mode beyond medium rota- 
tion rates in all sequences. Nevertheless, the ^/-mode fre- 
quencies still agree well quantitatively. Note however that 
avoided crossings of the ii"i-mode and the *pi-mode (as in 
sequences A and B) complicates a comparison. If the con- 
tinuous frequency curve of the ffi-mode is followed, the Hi- 
mode takes over the characteristics of the ^pi-mode at the 
avoided crossing. Thus, if continuous frequency curves are 
compared, a possible change of the mode labeling must be 
taken into account. Although the simulations in the Cowl- 
ing approximation reproduce the Hi, and ^pi-mode fre- 
quencies from the CFC simulations fairly close for most se- 
quences, which supports the idea of establishing empirical 
relations for predicting the correct frequencies from models 
evolved in the Cowling approximation, such relations should 
be used with caution (see discussion below). 

For sequence BU a comparison of models evolved us- 
ing a coupled evolution of hydrodynamics and spacetime, 
with no approximation for the gravitational field equations, 
against the Cowling approximation was presented by Font 
et al. (2002)^. Comparing these results with our new simu- 

^ Note that in that work the sequence already ended at our model 



Table 7. Relative difference in per cent in the frequencies of the 
-F, Hi, and ^pi-mode for the nonrotating model and the most 
rapidly rotating models of sequence A, AU, B, and BU in the 
Cowling approximation with respect to the CFC simulations. A 
-I- sign (— sign) indicates an overestimation (underestimation) of 
the CFC frequency in the Cowling approximation. 



Model 


F 


Hi 




^Pi 


AO/AUO/BO/BUO 


+86 


+15 


+16 


+10 


AlO 


+108 


-12 


-12 


+18 


AU5 


+77 


+5 


-9 


+22 


B9 


+81 


-1 


-3 


+18 


BU9 


+98 


+6 


-7 


+6 



lations using the CFC approximation as presented in Fig. 7 
shows excellent agreement for the frequencies of the _F-mode. 
This is another important confirmation of the validity of the 
results obtained with our code, considering in particular the 
differences in coordinate choice, grid resolution, evolution 
time, and assumption of CFC compared to the code em- 
ployed in Font et al. (2002)^ 

In the case of the //i-mode we also find rather good 
agreement, with increasing mismatch as the models ap- 
proach the mass-shedding limit. The oscillations in the Hi- 
mode frequency at large rotation rates observed in the fully 
relativistic simulations by Font et al. (2002) were interpreted 
as possible effects of avoided crossings. In contrast to this, in 
our simulations we find neither these oscillations nor signs 
of avoided crossing of the i/i-mode with other modes in se- 
quences AU and BU (see Fig. 7, and also Figs. 1 and 2). 
For sequences A and B, where there is clear evidence for 

BU7 with an axis ratio Vp/rc = 0.650, while several additional 
models with intermediate rotation rates were also evolved. 

A similar comparison for sequence AU yields equally good re- 
sults (Bernuzzi & De Pietri 2006). 
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Figure 7. Comparison of tlie frequencies for tlie F and _ffi-mode 
in CFC (black solid lines) and full general relativity (red dashed 
lines) for sequence BU. 



one in Font et al. (2002). Based on the work by Yoshida & 
Kojima (1997) and Yoshida & Eriguchi (2001), an empir- 
ical relation was derived using information from the com- 
pactness of a model, and assuming nearly linear scaling of 
the frequencies with increasing rotation rate. When applying 
this relation to our model AlO, we find a significant discrep- 
ancy between the predicted frequencies {fp ~ 0.882 kHz and 
fHi = 0.757 kHz) and the ones extracted from the actual nu- 
merical simulations (/f = 0.685 kHz and fni = 1.098 kHz), 
with relative differences of 29 and 31 per cent for the F and 
^/-mode frequencies, respectively. Although their assump- 
tion of nearly linear scaling of fp and fui with the rotation 
rate is confirmed by our results (see Figs. 3 to 6), this error 
is much larger than the predicted uncertainty of the relation 
of a few per cent. 

Owing to this, we refrain from establishing similar re- 
lations for other mode frequencies, even in cases where we 
also observe such nearly linear scaling of the frequency with 
rotation rate in CFC simulations and/or constant difference 
between the frequency of a specific mode in CFC and Cowl- 
ing. We rather suggest that except in special cases, (the 
computationally expensive) simulations of rotating stars in 
which the spacetime and the hydrodynamics are coupled 
cannot be replaced by a combination of simulations in the 
Cowling approximation and simple empirical relations. 



avoided crossing, this happens at larger rotation rates. We 
thus conclude that the oscillations visible in Fig. 7 (upper 
red dashed line) are an artifact which can be explained by 
the effectively lower grid resolution and shorter evolution 
time compared to the current simulations, as this increases 
the error in the numerical extraction of the mode frequency 
particularly for higher order modes. 

When comparing the frequencies for sequence BU from 
simulations with a dynamical spacetime to the ones in the 
Cowling approximation. Font et al. (2002) also observed that 
the frequencies for the -F-mode (and, to a lesser degree, also 
the ones for the //i-mode) show a similar dependence on 
the rotation rate in both cases. From the approximate con- 
stancy of the difference between the -F-mode frequency and 
the corresponding result in the Cowling approximation, they 
constructed an approximate empirical relation for calculat- 
ing the frequency of the _F-mode for arbitrary rotation rates. 
This relation thus depends on the F-mode frequency of the 
nonrotating model BUO obtained in a dynamical spacetime 
evolution and the variation of the F-mode frequencies with 
increasing rotation in the Cowling approximation. It yields 
the correct frequencies with an accuracy of better than 2 
per cent for the most rapidly rotating in their model series 
of sequence BU (our model BU7) . 

Unfortunately our results clearly suggest that it is not 
possible to straightforwardly apply this simple method to 
set up analogous relations for sequences A, AU, and B, as 
shown in Figs. 3 to 6. Only for sequence BU, the F-mode 
frequency exhibits a nearly constant difference between the 
CFC and Cowling simulations irrespective of the rotation 
rate. Moreover, the frequencies of all other modes displayed 
in these figures also do not fulfill such a relation, except 
maybe in the case of the ^pi-mode for sequences A and AU. 

Lacking results from simulations involving evolution of 
the spacetime, SAF proposed an approximate relation for 
the F and ^/-mode frequencies of sequence A similar to the 



5 EIGENFUNCTION RECYCLING 

As pointed out in Section 2.4, initial perturbations of the 
form given by Eqs. (11, 12) excite not only one desired spe- 
cific eigenmode, but additional oscillation modes as well. 
This is clearly visible in the power spectrum of model Al 
displayed in the upper panel of Fig. 8. While the funda- 
mental modes F and ^/ clearly dominate the spectrum, the 
overtones Hi and ^pi along with several other modes are also 
significantly excited by a perturbation using the I = and 
I — 2 trial eigenfunction of the initial velocity components 
Vr and ve, respectively. 

This consequence of using trial eigenfunctions can be 
avoided by performing an additional 'recycling' run. For this 
the actual two-dimensional eigenfunctions of the Vr and vg 
velocity components at a selected mode frequency are ex- 
tracted from the original simulation which was perturbed 
by a trial eigenfunction. These eigenfunctions are then ap- 
plied with appropriate amplitude as initial perturbation to 
a second simulation of the same stellar equilibrium model. 
In this case, other modes than the selected one are strongly 
suppressed, as shown for model Al in the lower panel of 
Fig. 8. Here eigenfunctions of the F and ffi-mode, and the 
^/ and ^pi-mode are extracted from a previous simulation 
perturbed with I — and I — 2 trial eigenfunctions, respec- 
tively. These are then used as initial perturbations for four 
recycling runs with relative amplitudes with respect to the 
amplitude of the F eigenfunction of 5.9, 3.8, and 7.0 for the 
H\, ^f, and ^pi eigenfunction, respectively, in order to arrive 
at approximately equal strengths of the dominant modes in 
the power spectrum. 

We note that the choice of using Vr and ve eigenfunc- 
tions as initial perturbations for doing the recycling simu- 
lation is not arbitrary. Naively, one would think to use all 
four evolved variables p, Vr, vg, and Vy, instead. However, 
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Figure 8. Fourier transform of the evolution of the radial profile 
of the rest-mass density p along the equatorial plane for model 
Al. For the two simulations in the upper panel an Z = and I = 2 
trial eigenfunction is used as perturbation, respectively, while the 
four recycling runs in the lower panel are excited with recycled 
eigenf unctions. The scaling of the ordinate is linear, and the power 
spectra in the upper panel are scaled to obtain the same strength 
of the _F-mode. 
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Figure 9. Radial profiles of the eigenfunctions for the Hi- 
mode (along the equatorial plane; upper panel) and for the 
^pi-mode (along 9 = 7r/4; lower panel) obtained from the evo- 
lution of model Al. For an initial perturbation using either a 
trial eigenfunction (black solid lines) or an extracted eigenfunc- 
tion (blue dashed lines), the shapes of the eigenfunctions agree 
well. The eigenfunctions are scaled to the same maximum height. 



for a given mode some of these variables are out of phase 
with respect to the others. Using the eigenfunctions of all 
variables as recycling perturbations simultaneously without 
taking into account the relative phase between them does 
not lead to the excitation of a single mode, but to the ex- 
citation of a sum of different modes (similar to choosing a 
trial eigenfunction). From the phase information contained 
in the complex FFT of the various variables, we determined 
that, at least for the modes we are interested in, the quanti- 
ties Vr and V0 have the same phase, while the other two are 
out of phase by 7r/2 with respect to and vg. 

While the suppression of undesired additional modes in 
the power spectrum significantly improves when eigenfunc- 
tion recycling is performed, the form of the eigenfunction of 
various metric and hydrodynamic quantities extracted from 
the recycling run is usually altered only negligibly as com- 
pared to the eigenfunction extracted from the original simu- 
lation with the trial perturbation. In Fig. 9 we present radial 
profiles of the rest-mass density eigenfunction p'^^ along the 
equatorial plane (upper panel) and of the ^-velocity eigen- 
function v'f along S = 7r/4 (lower panel). They are extracted 
from both the original simulations of model Al using I = 
(for p"^) and I = 2 (for vf) perturbations with trial eigen- 
functions, and also from the respective Hi and ^pi recy- 
cling runs. Only at the outer stellar boundary the shape of 
the eigenfunctions depends slightly on whether the model is 
perturbed by a trial or extracted eigenfunction, while in the 
bulk of the star the eigenfunctions are practically identical. 
We find similar results for the fundamental modes F and 
and for the eigenfunctions of other metric and hydrody- 
namic quantities in model Al and several other moderately 
rotating models. 



We can thus conclude that for such rotating stellar mod- 
els an initial perturbation with trial eigenfunctions is ade- 
quate to precisely obtain the mode frequency and to extract 
an accurate corresponding eigenfunction of the fundamental 
modes and their first overtones from the evolution. Addi- 
tionally, for such models a single recycling run suffices to 
efficiently suppress the excitation of all unwanted modes. 
However, if the peak of the investigated mode in the power 
spectrum is small and/or if several modes interact (see also 
the discussion of avoided crossings in Section 4.2), which is 
typically the case for higher order modes in rapidly rotating 
models, another recycling loop may be necessary to clearly 
determine the mode frequency and eigenfunctions, and to 
channel most of the initial perturbation energy into a single 
oscillation mode. 

Particularly in the case of the avoided crossing of the 
H\ and the *pi-mode in sequence A and B (see also Figs. 1 
and 2), we use a modified recycling strategy for an accurate 
mode analysis. Starting from a model where the mode fre- 
quency and eigenfunctions can still be clearly determined, 
the next model in the sequence is perturbed with the eigen- 
functions of the investigated mode extracted from the previ- 
ous model. This sequential recycling is a very helpful tool to 
resolve problems with ambiguous or unclear mode frequen- 
cies and character of the eigenfunction. 



6 NON-LINEAR PULSATIONS 

Although linear perturbations of rotating stars are assumed 
to have a vanishingly small amplitude, so that the back- 
ground equilibrium star is unaffected by a linear oscillation 
mode, in certain situations non-linear effects can become 
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Figure 10. Fourier transform of the evolution of the radial pro- 
file of the rest-mass density p along the equatorial plane for the 
nonrotating model (AO/AUO/BO/BUO) with an / = trial eigen- 
function used as initial perturbation. Additional to linear modes, 
several non-linear haj-monics can be identified. 

important. Our non-linoar evolution code allows us to inves- 
tigate such eflFects, of which we present several cases in the 
following. 

6.1 Non-linear harmonics 

The most basic non-linoar effect wo sec in our simulations 
is the appearance of non-linear harmonics of the linear pul- 
sation modes, a general property of non-linear systems (cf. 
Landau & Lifschitz 1976, § 28). To lowest order, these arise 
as linear sums and differences of different linear modes, in- 
cluding self-couplings. If the system has cigcnfrcqucncics LOi, 
the non-linearity of the equations excites modes at frequen- 
cies uii ± , with amplitudes proportional to the product 
of the amplitudes of the combining frequencies. We note 
that such non-linear harmonics were also recently noticed 
by Zanotti et al. (2005) in the numerical investigation of the 
dynamics of oscillating, relativistic, high-density tori around 
Kerr black holes. 

In Fig. 10 we present the Fourier PSD of the den- 
sity evolution of the nonrotating model (AO/AUO/BO/BUO), 
to which a finite radial I = initial perturbation of the 
form (11) was added. In addition to the main linear modes 
F and H\, one can observe several of their non- linear har- 
monics, such as the self-couplings 2 ■ F, 3 ■ F and the linear 
sums Hi — F and Hi — 2 ■ F. \n the PSD a large number 
of additional peaks can be seen, and essentially all of those 
should correspond to non-linear harmonics of the excited lin- 
ear modes. It is interesting to note that several non-linear 
harmonics arc present that have frequencies much smaller 
than the fundamental radial mode F (which, in the linear 
approximation, possesses the lowest frequency). These can 
actually fall into the frequency range of the inertial modes 
for rotating models. Thus, in rotating models further non- 
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Figure 11. Same as Fig. 10, but with an ( = 2 trial eigenfunction 
used as initial perturbation. 

linear interactions between radial modes and inertial modes 
can be expected (see also the discussion in Section 6.2). 

Non-linear harmonics can occur not only due to cou- 
plings between modes of the same index I (such as the I = 
modes discussed above), but also duo to couplings between 
modes of different I. Fig. 11 shows several identified harmon- 
ics that represent linear sums and differences between linear 
modes for the nonrotating model, when an Z = 2 perturba- 
tion was added to the initial data of the nonrotating model. 
Due to the approximate nature of the chosen trial eigenfunc- 
tion, I = modes are also excited in addition to the main 
I — 2 modes. The presence of both I = and 1 = 2 modes 
then leads to the appearance of several non-lineax harmon- 
ics, which include cases like '^pi — F, Hi — etc. It is thus 
clear that even though in the linear approximation modes 
of different / are orthogonal to each other (in a nonrotating 
perfect fluid star), non- linear effects couple all linear modes 
with different I that are present in a non- linear simulation)^. 

In rotating models, where there can be a perplexing 
alternation of various linear modes and non-linear harmon- 
ics in a PSD plot, one can easily distinguish the non-linear 
harmonics from the linear modes by comparing PSDs pro- 
duced from simulations with different initial perturbation 
amplitudes. As expected, the linear modes scale almost lin- 
early, while the non-linear harmonics scale as the product of 
the amplitudes of the modes (or of the modes and harmon- 
ics) from which they are produced. Such a case is shown in 
Fig. 12, where model Al from sequence A is used in which 
the linear modes F, Hi. '^f, and ^pi are all excited at roughly 
the same strength using the eigenfunction recycling tech- 
nique for the initial perturbation as described in Section 5. 
Comparing two simulations of this model that differ by a 
common factor of 4 in the initial perturbation amplitude. 



^ Notice, however, that for relativistic stars there exists no proof 
yet on the completeness of quasi-normal modes. 
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Figure 12. Fourier transform of the evolution of the radial profile 
of the 0- velocity Dg along 6 = 7r/4 for model Al with a combined 
-F, Hi, and ■^pi-mode eigenfunction as initial perturbation. 
Changing the perturbation amplitude a reveals the non-linear 
scaling of the coupled modes as opposed to the linear scaling 
of the linear modes. 

one can clearly notice that while the amplitudes of the lin- 
ear modes scale nearly linearly, the amplitudes of the various 
non-linear harmonics scale non-linearly (with many of them 
scaling roughly quadratically) . 

6.2 Non-linear 3-mode couplings 

The presence of non-linear harmonics opens the possibility 
of 3-mode couplings when the star is rotating. The reason 
such 3-mode couplings can take place is the fact that the ef- 
fect of rotation on the different modes varies. As was already 
discussed in Section 4.1, higher-order modes are typically af- 
fected stronger by rotation than lower-order modes, which 
results in avoided crossings between mode sequences. Rota- 
tion also influences the frequencies of the various non-linear 
harmonics to a different degree. Thus, at certain rotation 
rates, a non-linear harmonic can have the same frequency as 
a linear mode. Examples of such cases are demonstrated in 
Fig. 13, which shows the frequencies of several linear modes 
and of the two non- linear harmonics ^pi — and Hi — F as a 
function of T/j W| . Apparently the ^pi —^f harmonic is cross- 
ing the frequency of the fundamental quasi-radial F-mode 
at about T/\W\ ~ 0.1, while the Hi — F harmonic is crossing 
the frequency of the ^pi-mode at about ~ 0.12. We 

also note that at any rotation rate one can expect several 
non-linear harmonics to coincide in frequency with some of 
the infinitely many inertial modes contained in the inertial 
mode range. 

At such crossings, the coinciding of the frequencies 
could potentially lead to resonance effects and even to para- 
metric instabilities. In such a way, significant energy from 
one mode could be transfered to other modes. The most in- 
teresting case would be if pulsational energy from the quasi- 



Figure 13. Frequencies of various linear modes and the two non- 
linear harmonics Hi — F and ^pi — -^f for sequence A. At specific 
rotation rates the frequencies of linear and non-linear modes co- 
incide, which can lead to non-linear 3-mode coupling. 



radial mode, which weakly radiates gravitational waves, 
could be transfered non-linearly to stronger radiating non- 
radial modes. Since during core bounce a significant amount 
of kinetic energy is stored in the radial modes of pulsation, 
the transfer of even a small percentage of this kinetic energy 
to a nonradial mode could result in the emission of strong 
gravitational waves. This scenario has first been suggested 
by Clark (1979) and Eardley (1983), who investigated the 
parametric instabilities that could take place, using a New- 
tonian, slowly rotating collapse model (for recent related 
work for nonrotating or slowly rotating relativistic stars, see 
Passamonti et al. 2005). 

In nonrotating or slowly rotating models, such a para- 
metric instability can only take place under special condi- 
tions that would allow the two modes to be in resonance. 
Here we find that in rapidly rotating models rotational shift- 
ing of the frequency of different modes broadens the range 
of parameters for which interesting resonances could take 
place. We particularly notice that the quasi-radial mode will 
be in resonance with some inertial mode(s) for all rotation 
rates above a critical value. It is thus interesting to further 
study the possible energy transfer between different modes 
excited after, e.g., a core collapse or an accretion-induced 
collapse event, either on secular time-scales or as a para- 
metric instability. 

Here we only observe the necessary conditions for non- 
linear 3-mode couplings to take place*. Whether such cou- 
plings will indeed lead to strong parametric resonances and 
enhanced gravitational wave emission remains to be investi- 
gated through much more detailed studies. 



■* We include the case of self couplings, when discussing 3-mode 
couplings. 
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Figure 14. Distribution of the logarithm of the rest-mass den- 
sity, logp (g cm""^), in the meridional plane for model BU9 with 
ideal fluid EOS at t = 4.4 ms in the Cowling approximation. 
The ordinate coincides with the rotation axis. The single shocks 
traveling through the atmosphere and creating a matter envelope 
are clearly visible. The black color coding corresponds to the low 
density artificial atmosphere where the rest-mass density is actu- 
ally much smaller than 10^ g cm""^. Note that only the innermost 
80 km of the computational grid are shown. 



6.3 Mass-shedding-induced damping 

Another striking example of a non-linear effect is the mass- 
shedding-induced damping of oscillations in stellar models 
which rotate at or near the mass-shedding limit. This new 
damping mechanism was first observed and discussed by 
SAF, and is especially important for pulsations in uniformly 
rotating stars. In SAF, the mass-shedding-induced damping 
was demonstrated for the uniformly rotating models BUS 
and BU9 using fixed-spacetime evolutions, i.e. the Cowling 
approximation. 

As explained in SAF, the damping mechanism works 
as follows. As the star approaches the mass-shedding limit, 
the effective gravity near the equatorial surface diminishes, 
exactly vanishing at the mass-shedding limit. A small radial 
pulsation then suffices to cause mass-shedding after each os- 
cillation period. As a result, a low-density envelope is created 
around the star. This envelope is initially concentrated in the 
regions close to the stellar equator, but with each oscillation 
period more and more mass is shed in the form of shock 
waves, and the envelope expands outwards and away from 
the equatorial plane. Since in rotating stars every pulsation 
mode also has a radial velocity component, the damping af- 
fects all modes. In SAF, it was found that the damping in 
the Cowling approximation can be rather strong. Here we 
investigate the same damping effect in the CFC approxima- 
tion and compare the two cases. 

We first repeat the study of the mass-shedding-induced 
damping for the radially perturbed model BU9, as presented 



M 40.0 




40.0 
r [km] 



Figure 15. Same as Fig. 15, but for a spacetime coupled to the 
hydrodynamics (CFC simulation). Note that in contrast to the 
Cowling approximation, a much smaller angular part of the atmo- 
sphere is filled with a matter envelope, which also has significantly 
also lower density. 



in SAF, in the Cowling approximation. Fig. 14 shows the 
distribution of the rest-mass density p in the meridional 
plane using the (nonisentropic) ideal fluid EOS (10). At time 
t = 4.4 ms several oscillations have already occurred and the 
high-entropy envelope has acquired a near-equilibrium state, 
extending into almost the entire computational grid (whose 
outer boundary has been set to 5 stellar equatorial radii with 
120 additional logarithmically spaced radial grid points)^. 
Near the equatorial plane the envelope has a rest-mass den- 
sity of ~ 10~® to 10~*pc, through which further consecutive 
shocks propagate. Only at angles < 20° (as measured from 
the rotation axis) does the high-entropy envelope not com- 
pletely form. 

We then study the same perturbed model (with a sim- 
ilar initial effective perturbation amplitude) but also evolv- 
ing the spacetime, in the CFC approximation. Compared 
to the Cowling approximation, the mass-shedding behavior 
and the properties of the matter envelope are now signifi- 
cantly different. In Fig. 15, one can see (at the same time 
t — 4.4 ms as in Fig. 14) that the high-entropy envelope 
outside the star is confined to within roughly 45° with re- 
spect to the equatorial plane, filled with matter of one to 
two orders of magnitude lower density than in the Cowling 
approximation. In addition, the consecutive shocks barely 
reach the outer grid boundary. 

Fig. 16 shows the corresponding radial rest-mass density 
profiles along the equatorial plane for both the Cowling and 
CFC approximations. The envelope in the CFC approxima- 



^ Note that the equatorial stellar radius is at re ~ 17 km. In 
Table 1 dimensionless units are used for , resulting in a different 
numerical value of the same physical location. 
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Figure 16. Radial profiles of the rest-mass density p along the 
equatorial plane for model BU9 with ideal fluid EOS at t = 4.4 ms 
in CFG (black solid line) and in the Cowling approximation (red 
dashed line). The single shocks due to mass-shedding are clearly 
visible. The dotted line marks the equatorial stellar radius. 



tion has a lower density than in the Cowling approximation, 
with the difference consistently increasing with distance and 
becoming 2 orders of magnitude at the outer grid bound- 
ary. In this plot the single shocks generated by the stellar 
oscillations, which first create the envelope and later prop- 
agate through it, are clearly noticeable. For model BU9 the 
frequency of the F-mode oscillations, which are mainly re- 
sponsible for the mass-shedding in the above simulations is 
about twice as high in the Cowling than in the CFC approx- 
imation (2.313 kHz vs. 1.169 kHz). This is, however, not re- 
flected in the shock pattern of the radial density profile in 
Fig. 16, where naively twice as many shocks between Vc and 
some radius r > could be expected. This is a consequence 
of the unequal profiles of the shock propagation velocity in 
the envelope due to the different density profiles, as well as 
interactions of shocks traveling at different velocities due to 
unequal shock strengths. The latter effect can also lead to 
reverse shocks, which further complicate the picture. Note 
that in the density profile from the simulation in the Cowling 
approximation in Fig. 16, the significant loss of the matter 
in the stellar regions close to re by mass-shedding is evident, 
whereas in the CFC simulation the initial density profile in- 
side the stellar boundary is nearly preserved. 

In an equilibrium model rotating at the mass-shedding 
limit, effective gravity at the equatorial surface vanishes due 
to a delicate balance between the pressure gradient, the 
gravitational force, and the centrifugal force. In the Cowling 
approximation the spacetime is held fixed during the entire 
simulation, and thus the gravitational force cannot react to 
local over- or underdensity caused by the pulsations. Evi- 
dently, in this case matter can easily be ejected from the 
stellar equatorial surface at each pulsation. In contrast to 
this, when the spacetime is coupled to matter it responds to 
changes in the hydrodynamic evolution. The most strongly 
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Figure 17. Radial profiles of the _F-mode eigenfunction of the 
rest- mass density (black solid line), the conformal factor ijf^ 
(blue dashed line), and its negative radial derivative —dr<f>''^ (red 
dotted line) along the equatorial plane for model BU9 with ideal 
fluid EOS in CFC. The vertical dotted lines mark the locations 
^np, ^n^^ ^ndr-il> °f eigenfunction nodes and the equatorial 
stellar radius re. The vertical scaling of the eigenfunctions is ar- 
bitrary, but their sign is not. 

excited pulsation mode in the mass-shedding simulations of 
model BU9 is the F-mode, whose rest-mass density eigen- 
function p'^^ features a node at rn p — 7 km in the equatorial 
plane®, located at less than half the stellar equatorial ra- 
dius (see Fig. 17). However, not everywhere in the star an 
over- or underdensity is reflected by an according localised 
increase or decrease of gravity, as the local gravitational pull 
is created by the global density distribution. Consequently, 
the nodes of the eigenfunctions of both the conformal factor 
0°' and its negative radial derivative —dr4>^^ (representing 
the perturbation of the gravitational force) lie much farther 
out at rnd^tj, — 14 km > rn0 — 11 km > rnp. In a region 
rnp < r < r^^ the perturbations of density {p"^) and gravi- 
tational force (0"^, —drcj)"^) have opposite sign, and thus at 
some point during a pulsation cycle the gravitational force 
there decreases although the local density increases. Never- 
theless, close to the stellar equatorial surface both cf>''^ and 
—drijf^ are again in phase with p"^, and remain greater than 
zero (albeit rather small) at the surface. Thus, whenever the 
local density there increases in a pulsation, this small local 
net gain in the gravitational pull suffices to efficiently atten- 
uate mass-shedding, as observed in the CFC simulations. 

When matter is shed from the star, kinetic energy is 
carried away to the expense of the pulsational energy. Con- 
sequently, the pulsations of the star are gradually damped, 
depending on the intensity of mass-shedding. Fig. 18 shows 
the time evolution of the central rest-mass density for 
model BU9 obtained from the CFC simulation and the one 



° The occurence of at least one node in the radial eigenfunction 
of a quasi-radial mode is a consequence of mass conservation. 
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Figure 18. Time evolution of the central rest-mass density pc 
for model BU9 with ideal fluid EOS at t = 4.4 ms in CFC (black 
solid line) and in the Cowling approximation (red dashed line). 
In the Cowling approximation, significant mass-shedding-induced 
damping can be seen, whereas the decline of the oscillation am- 
plitude in CFC can be mostly attributed to mode-bcating. 



in the Cowling approximation. It is evident that the damp- 
ing of the pulsations is very strong in the latter case, leading 
to an approximate damping time-scale of 10 ms. In contrast, 
the insignificant mass-shedding in the CFC simulation re- 
sults in only small damping on a much longer time-scale. 
Actually, a long-term observation of the time evolution of 
Pc reveals that most of the decline of the oscillation ampli- 
tude in the time window of Fig. 18 can be explained by a 
mode-beating effect, which superimposes the genuine damp- 
ing. Even with high initial perturbation amplitudes we are 
not able to observe significant unambiguous mass-shedding- 
induced damping during typical evolution times of 20 ms. 
Simulations with lower initial perturbation amplitudes or of 
models which rotate not so close to the mass-shedding limit 
exhibit less mass-shedding compared to the above model, 
which is in accordance with the findings presented by SAP. 

Even though in the CFC simulation the mass-shedding- 
induced damping of pulsations is not as strong as in the 
Cowling approximation, it could still have significant impli- 
cations for unstable modes that grow on secular time-scales 
due to the gravitational-radiation driven CPS instability. In 
particular, the I = m = 2 /-mode becomes unstable only 
near the mass-shedding limit in uniformly rotating stars 
(see Stergioulas & Friedman 1998; Morsink, Stergioulas & 
Blattnig 1999). A detailed investigation is required to deter- 
mine whether the rate of mass-shedding-induced damping 
is shorter than the secular growth rate of the instability at 
amplitudes smaller than of order unity. 

In Section 2.4 we have demonstrated that for exciting 
-F-mode oscillations, an initial velocity perturbation of the 
form (11) is more appropriate than the density perturbation 
of Eq. (14) if the spacetime is coupled to the hydrodynamics. 
However, we have performed mass-shedding test simulations 



of model BU9 in the Cowling approximation, which show 
that in this case a velocity perturbation with typical ampli- 
tudes either results in a strong negative or positive drift in 
the time evolution of the central rest-mass density, depend- 
ing on the sign of the initial perturbation amplitude. Avoid- 
ing this by significantly reducing the initial perturbation am- 
plitude in turn effectively suppresses the mass-shedding and 
the resulting damping of the stellar pulsations. Following 
SAF, we therefore use the density perturbation (14) for the 
simulations of the mass-shedding model BU9 in both CFC 
and the Cowling approximation presented here. To compen- 
sate for the amplification of the oscillation amplitude dur- 
ing evolution using equal initial perturbation amplitudes for 
a coupled spacetime as compared to the Cowling approxi- 
mation (see Section 2.4), we use different initial perturba- 
tion amplitudes in CFC (a = 0.01) than in the Cowling ap- 
proximation (a = 0.05). This yields an approximately equal 
height of the first oscillation peak during the evolution (see 
Fig. 18) , so that we can consider the initial conditions in the 
two simulations to be similar. 



7 GRAVITATIONAL WAVES 

In the numerical simulations of the models presented above, 
we extract the gravitational waves emitted by the pulsa- 
tions initiated in each model. For this we use the New- 
tonian quadrupole formula in its time-integrated form (as 
described in detail in Dimmelmeier, Font & Miiller 2002b), 
which yields the quadrupole wave amplitude Afo as the low- 
est order term in a multipole expansion of the radiation field 
into pure-spin tensor harmonics (Thorne 1980). 



7.1 Gravitational wave power spectrum 

Fig. 19 shows the gravitational wave power spectrum (i.e. 
the PSD of the time evolution of Afo ) for a series of simula- 
tions of model Al (i.e. the most slowly rotating model of se- 
quence A) in which the F, Hi, ^f, and ^pi-modes are individ- 
ually excited using the eigenfunction recycling technique. As 
already stated in Section 5 the initial excitation amplitude 
a for each individual mode was chosen such that the PSDs 
of the integrated density variations in the equatorial plane 
reach similar amplitudes (see also lower panel of Fig. 8). 
However, this choice for a results in a different relative os- 
cillation amplitude Spc of the central density during the evo- 
lution (with 2Spc being the density variation measured top 
to bottom of the oscillation) for each of the four simulations. 
In Fig. 19 we thus scale the original gravitational wave signal 
such that it corresponds to 5pc = 1 per cent. This is possi- 
ble as we find that the maximum amplitude [A^q |niax of the 
quadrupole signal is approximately proportional to Spc (and 
in turn also to a) to fair accuracy for the four linear stellar 
pulsation modes investigated here. 

As expected, the quadrupolar ^/-mode is the strongest 
emitter of gravitational waves, while the other three modes 
are still within roughly one or two orders of magnitude of 
this dominant gravitational wave emitting mode. Note that 
an excitation of the ^/-mode results in a nearly monochro- 
matic gravitational wave power spectrum, indicating that 
the recycled eigenfunctions used for exciting this mode are 
of high accuracy. The PSD for the individual excitations of 
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Figure 19. Gravitational wave power spectrum for model Al 
excited with recycled eigenfunctions of the F, Hi, and ^p\- 
mode, respectively. In the F, Hi, and ■^pi-mode recycling runs 
a significant ^/-mode component is also present in the gravita- 
tional wave signal. The thick lines indicate an artifical fall-off 
with quadratic frequency dependence assumed to get rid of these 
spurious contributions. Units on the i/-axis are arbitrary. 



the F, Hi, and ^pi-mode are also dominated by the ap- 
propriate frequency of the respective mode, but exhibit also 
additional contributions of other excited modes. These con- 
tributions are usually small (at least one order of magnitude 
smaller than the dominant mode), but as the quadrupole ^J- 
mode strongly radiates gravitational waves, a small leaking 
of energy into that mode from the main excitation mode is 
sufficient to exhibit a significant contribution at that specific 
frequency. This is particularly apparent in the case of the 
//i-mode, which radiates gravitational waves only through 
its rotationally-induced I ^ 2 terms. 

In order to suppress the unwanted '^/-mode contribution 
in the computation of the characteristic signal frequency and 
amplitude in the following Section 7.2, for the gravitational 
wave power spectrum of the F, Hi, and ^pi-mode simu- 
lations we have assumed a fall-off behavior with quadratic 
frequency dependence in the region of the spectrum where 
the spurious ^/-mode is located (indicated by thick lines in 
Fig. 19). 

We also note that for the excitations of the overtones 
Hi and ^pi with recycled eigenfunction, some inertial modes 
(visible at the low-frequency end in the PSD of Fig. 19) also 
weakly contribute to the emission of gravitational waves. 



7.2 Detectability of gravitational waves 

In order to assess the prospects for detectability of the grav- 
itational wave signal from our models of pulsating neutron 
stars by current and planned interferometer detectors, we 
calculate the characteristic dimensionless amplitude he of 
the signal from the quadrupole wave amplitude A20 as de- 



Table 8. Characteristic frequency fc and characteristic ampli- 
tude he for a gravitational wave signal emitted by the most slowly 
rotating model Al excited with recycled eigenfunctions of the _F, 
-ffi, and ^pi-mode, respectively, as prospectively seen by var- 
ious interferometer detectors. For comparison, the frequency / of 
the respective stellar oscillation mode is also given. All frequen- 



cies are 


in kHz, 


and he 


is given 


in units of 10 


■^■^ for a source 


duti T'gw 


= 10 kpc oscillating with Spc = 1 per 


cent. The signal 


duration time is 


assumed to be t^. 


rw = 20 ms. 








VIRGO 


LIGO I 


adv. LIGO 


Mode 


/ 


fc 


he 


fc he 


fc he 


F 


1.400 


1.370 


19.3 


1.338 18.9 


1.236 17.4 


Hi 


3.816 


3.783 


12.1 


3.792 12.1 


3.791 12.1 




1.577 


1.524 


463.3 


1.467 467.5 


1.292 393.8 




3.580 


3.545 


76.9 


3.545 76.9 


3.544 76.9 



scribed in Zanotti et al. (2005). We perform a Fourier trans- 
form of the transverse traceless gravitational wave signal, 

/oo 
e^^'f'h'^{t)dt, (15) 

where 

^TT^l /15^ (16) 

8 V TT Tgw 

in the equatorial plane of the neutron star (assuming optimal 
detection geometry of the interferometer) with rgw being the 
distance from the emitting source to the detector. To obtain 
the detector dependent characteristic signal frequency 

and characteristic signal amphtude 

he ^[3 ^{\h\'}fdfj, (18) 

the power spectral density Sh of the detector is needed 
(with She = Shifc))- We approximate the average {\h\^) 
over randomly distributed angles by \h\^. From Eqs. (17, 
18) the signal-to-noise ratio can be computed as S/N — 
/ic/[/irms(/c)], where /irms(/c) = ^/ feSh{fe) IS the valuc of 
the rms strain noise (i.e. the theoretical sensitivity window) 
for the detector at the characteristic frequency. 

In Table 8 we give the values of fe and he for the grav- 
itational wave signal emitted by the most slowly rotating 
model Al excited with recycled eigenfunctions, computed 
for the theoretical power spectral density Sh of the detectors 
VIRGO, LIGO I and advanced LIGO located at a distance 
of r-gw = 10 kpc from the source. The strong emission of the 
gravitational waves by the quadrupolar ^/-mode and its first 
overtone ^pi is reflected by the large wave amplitudes he of 
~ 400 X 10~^^ and ~ 80 x 10~^^, respectively, compared to 
only ~ 20 X 10"^^ for the F-mode and ~ 10 x 10~^^ for the 
Hi-moAe. 

Note that we integrate /i^^ over a finite duration time 
of tgw = 20 ms in the FFT and scale 5pe to 1 per cent for 
each individual excited mode. To compute he in a straight- 
forward way from the values given in Table 8 for an arbitrary 
source distance, signal duration, and central oscillation am- 
plitude of the star, the following simple scaling laws can be 
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Figure 20. Distribution of the characteristic gravitational wave 
signal amplitude he in the frequency space for model Al excited 
with recycled eigenfunctions of the F, Hi , ^f, and ^pi-mode, re- 
spectively. The crosses represent he for sl source a,t ^'gw — 10 kpc, 
oscillating with 5pc = 5 per cent. The signal duration time is 
assumed to be tgw = 20 ms. The sensitivities hrms of the interfer- 
ometer detectors VIRGO (black solid line), LIGO I (red dashed 
line), and advanced LIGO (blue dotted line) are also shown. Par- 
ticularly for the F and ^/-mode signal, the value of he and fc 
is slightly dependent on the sensitivity of the respective inter- 
ferometer detector, as indicated by the color coding. The arrows 
indicate the location of he for 5pc = 1 per cent. 



used: he 



and he 



while he cc \A 



E2| 

20 I max 



for 



a given power spectral density Sh and thus approximately 
proportional to Spe (see discussion in Section 7.1). 

As we artificially remove the spurious contribution of 
the quadrupole ^/-mode to the gravitational wave power 
spectrum of the F, Hi, and ^pi-mode recycling runs (see 
Section 7.1), the small discrepancy between the values for 
characteristic gravitational wave signal frequency fe and the 
stellar oscillation mode frequency / in Table 8 is predomi- 
nantly caused by the dependence of fc on the power spectral 
density Sh of the specific detector. Accordingly, as he 7^ h'^'^ , 
the values of characteristic amplitudes he are slightly dif- 
ferent for each detector. The effect of Sh on fe and he is 
strongest for the gravitational wave signal emitted by the 
relatively low frequency modes F and 

In Fig. 20 we show the location of the characteristic 
gravitational wave amplitude he for the F, Hi, ^f, and ^pi- 
mode signal of model Al in relation to the sensitivity win- 
dows of VIRGO, LIGO I, and advanced LIGO. We again 
assume a distance to the source of rgw = 10 kpc, i.e. within 
our Galaxy, and (in contrast to Table 8) an oscillation ampli- 
tude Spe = 5 per cent. In such a scenario the gravitational 
wave signal from a neutron star pulsating in the ^/-mode 
with the adopted oscillation amplitude would be clearly de- 
tectable even by current interferometer observatories with a 
signal-to-noise ratio of around 50 for VIRGO and LIGO I 
(and 300 for the planned advanced LIGO detector). On the 
other hand, the F and ^pi-mode signals lie only marginally 



above the sensitivity threshold of current detectors, corre- 
sponding to a signal-to-noise ratio of at most ~ 2. The 
wave signal from the _ffi-mode, in turn, could possibly be 
detectable only by the proposed advanced LIGO detector 
for the selected pulsation amplitude and duration. 

If Spe is reduced to 1 per cent, as exemplified by the ar- 
rows in Fig. 20, the F and '^pi-mode signals clearly fall out 
of the sensitivity window of VIRGO and LIGO I, but would 
be detectable by advanced LIGO, while the _ffi-mode signal 
is now evidently out of the range even for the latter. Only 
the signal from the ^/-mode still exhibits a signal-to-noise 
ratio of around 10 for the current detectors. It is impor- 
tant to note that using the scaling properties of the signal 
discussed above, the data presented in Fig. 20 can be used 
to determine the necessary minimum stellar oscillation am- 
plitude Spe, maximum source distance rgw, and minimum 
signal duration time fgw required for a prospective detec- 
tion (set by a specific signal-to-noise ratio above 1) of the 
gravitational wave emitted by model Al. 

We also emphasize that the frequencies of the investi- 
gated linear pulsation modes F, Hi, ^f, and ^pi are located 
at the rapidly ascending high frequency slope of the detector 
sensitivity curves. Consequently, for the fixed rest mass se- 
quences A and AU the decrease of the frequencies for these 
modes with increasing rotation can enhance the prospects 
for detection by shifting the signals towards the more sensi- 
tive domain of the detector sensitivity window. For the fixed 
central density sequences B and BU, however, the opposite 
dependence of the Hi and particularly of the ^/-mode fre- 
quencies on rotation reverses this effect for these two modes. 

Further increase of the signal-to-noise ratio of the grav- 
itational wave signals presented here could be achieved with 
the advanced LIGO detector, using its planned ability to 
be tuned at higher frequencies in a narrow-banding mode. 
The chances for a successful detection of neutron star oscil- 
lation modes will also rise dramatically by the simultaneous 
operation of high-frequency-band detectors such as the pro- 
posed dual sphere detectors (Cerdonio et al. 2001), which 
are designed to have an improved sensitivity at frequencies 
of a few kHz and would thus be particularly suitable for 
detecting the first overtones Hi and ^pi. 



7.3 Gravitational wave asteroseismology 

It is clear from the discussion of our results so far, that pul- 
sating rotating neutron stars are gravitational wave sources 
that depend on several parameters (EOS, mass, angular mo- 
mentum, differential rotation law, initial perturbation am- 
plitude, damping mechanisms, etc.). All these parameters 
may have different effects on the oscillation spectrum of the 
star and, therefore, the successful extraction of the physical 
characteristics of the source from the gravitational wave sig- 
nal will be difficult to achieve. It is important to isolate each 
effect on the gravitational waveform in order to find general 
trends. 

Empirical formulas that can be used for gravitational 
wave asteroseismology (with the aim of pin-pointing the cor- 
rect high-density EOS through gravitational wave observa- 
tions) were constructed in recent years (Andersson & Kokko- 
tas 1998; Kokkotas, Apostolatos & Andersson 2001; Benhar, 
Ferrari & Gualtieri 2004), assuming nonrotating cold equi- 
librium models. Thus, these formulas contain only the effect 



© 2005 RAS, MNRAS 000, 1-24 



22 H. Dimmelmeier, N. Stergioulas and J. A. Font 



of EOS and mass. More recently, finite temperature effects 
and slow-rotation effects on the mode frequencies were com- 
puted (Ferrari, Miniutti & Pons 2003a,b; Ferrari et al. 2004), 
but these effects have not yet been incorporated in the em- 
pirical formulas for asteroseismology. In our work we focus 
on the effects of rotation, using a particular cold EOS as 
an example. A complete treatment will, of course, require 
using a large sample of realistic EOSs, but such an attempt 
is beyond the scope of this work. Using only the frequen- 
cies extracted from the simulations presented here, we can 
nevertheless identify several interesting features that indi- 
cate qualitatively the influence of rotation on the empirical 
formulas used in gravitational wave asteroseismology. 

First, the frequencies of most modes show a simple de- 
pendence on rotation so that the empirical formulas ob- 
tained for nonrotating stars could be extended to rotation by 
multiplying with a factor containing only one or two powers 
of a rotational parameter such as r/|VK|. The exception is 
the Jfi-mode, which will require special treatment due to the 
avoided crossing at moderate values of T/|iy|. Second, the 
frequencies of most modes, such as the F, and -ffi-mode, 
become nearly independent of the degree of differential ro- 
tation when parametrized as a function of r/|W|. However, 
the frequency of the ^pi-mode is strongly affected by the 
degree of differential rotation. Thus, am additional parame- 
ter, measuring the strength of differential rotation, could be 
built into an empirical formula for this specific mode. 

Ideally, if the four modes F, Hi, ^f, and ^pi were de- 
tected, then a set of empirical formulas for their frequencies, 
constructed in the way described above, would allow the si- 
multaneous extraction of the mass, radius, T/jWj, and de- 
gree of differential rotation. Of course, this ideal situation 
may become more complicated by the inclusion of finite tem- 
perature effects, which would add at least one more param- 
eter to the system. 



8 SUMMARY AND OUTLOOK 

Using the axisymmetric general relativistic hydrodynamics 
code CoCoNuT we have investigated pulsations of uni- 
formly and differentially rotating neutron star models. We 
have compared our numerical simulations, in which the 
spacetime dynamics is coupled to the evolution of the fluid, 
to previous results performed under the assumption of a 
fixed spacetime (Cowling approximation). In the present 
work we have used the so-called conformal flatness condi- 
tion (CFC) for the spatial three-metric, this being an ex- 
cellent and well-tested approximation of the exact space- 
time in the regime studied here. The coupled system of (hy- 
perbolic) hydrodynamics and (elliptic) metric equations has 
been solved using the novel approach of combining Riemann- 
solver-based HSRC methods for the fluid evolution and spec- 
tral methods for the computation of the spacetime met- 
ric (Dimmelmeier et al. 2005). 

As equilibrium initial data we have constructed four 
sequences of relativistic polytropes with parametrized rota- 
tion, which have then been perturbed by small amplitude 
Z = 0, 2, or 4 trial eigenfunctions to excite pulsation modes. 
By analyzing the time evolutions of various hydrodynamic 
and metric quantities using Fourier transforms along ra- 
dial profiles or on the entire spatial grid, we have obtained 



the pulsation frequencies for the fundamental quasi-radial 
{I = 0) F-mode, the fundamental quadrupole [l — 2) '^f- 
mode, and their respective overtones, the Hi and '^pi-mode, 
as well as three inertial modes labelled i_2, ti, and 12. Ad- 
ditionally, for two differentially rotating sequences we have 
obtained the frequencies for the and "^pi-mode. We have 
found that for these two sequences the "'pi -mode engages 
in an avoided crossing with the f/i-mode. This is a conse- 
quence of the different influence of rotation on the frequency 
of these two modes, which brings their frequencies closer 
with increasing rotation rate. At the avoided crossing, an 
exchange of the character of eigenfunctions of the Hi and 
*j3i-mode has been observed. Although linear perturbation 
theory predicts the existence of an infinite number of iner- 
tial modes in a finite frequency range, in our simulations 
we have detected the predominant excitation of only a few 
specific inertial modes, which are excited as by-products of 
the excitation of other modes. For the differentially rotating 
sequence of fixed rest mass, the inertial mode frequencies 
reach a maximum value before the mass-shedding limit is 
reached and then decrease again. 

In order to suppress the simultaneous excitation of more 
than one mode, which is typically the CEise when trial eigen- 
functions are used as initial perturbations, we have employed 
a new technique, which we call eigenfunction recycling. For 
this, we first extract the two-dimensional eigenfunction of 
velocity components at the frequency of the mode under 
consideration and then apply it as an initial perturbation 
in a second simulation of the original neutron star model. 
This selective excitation of modes works very well even for 
overtones. 

When comparing the frequency curves of the F-mode 
for all four investigated sequences with previous results in 
the Cowling approximation, we have observed that the lat- 
ter typically overestimate the correct frequency by about 
a factor of 2 (corresponding to an absolute difference of 
~ 1 kHz), which is consistent with similar findings for non- 
rotating stars by Yoshida & Kojima (1997). Moreover, we 
have found no evidence for the F-mode splitting found by 
SAF to occur in the Cowling approximation. Thus, we con- 
clude that this effect is an artifact which can be attributed 
to the Cowling approximation, as was suggested in SAF. On 
the other hand, for the H, and ^pi-mode frequencies, 
much closer agreement with the Cowling results of SAF has 
been found. For these modes the relative difference is usually 
less than ~ 20 per cent. Concerning empirical formulas con- 
structed by Font et al. (2002) and SAF using results in the 
Cowling approximation in order to predict mode frequen- 
cies for rotating sequencies, we have found that the chang- 
ing characteristics of equilibrium models along various se- 
quences significantly aifect the relative differences between 
the actual frequencies and those obtained in the Cowling ap- 
proximation. As a result, such empirical relations can only 
be of limited validity, while still representing a significant 
improvement compared to the Cowling approximation. 

Even though we have used only small-amplitude initial 
perturbations, the use of a non-linear evolution code has 
allowed us to investigate various non-linear effects in pul- 
sating neutron star models. We have observed non-linear 
harmonics of the linear pulsation modes of our stellar mod- 
els, arising as linear sums and differences of various linear 
modes, including self-couplings (for similar findings in high- 
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density accretion tori around Kerr black holes, see Zanotti 
et al. 2005). For rotating models such non-linear harmonics 
actually fall into the frequency range of the inertial modes 
and could interact with linear modes in non-linear 3-mode 
couplings in the form of a resonance or a parametric in- 
stability. It would be interesting to investigate under which 
conditions strong mode couplings or instabilities can occur. 

In SAF a new non-linear damping mechanism for non- 
linear pulsations was found to operate in uniformly rotating 
models that axe near the mass shedding limit. The damping 
of pulsations is due to mass shedding, when fluid elements 
near the equator become unbound due to the radial com- 
ponent of the velocity perturbation. This mass-shedding- 
induccd damping could have severe consequences for, e.g., 
the /-mode gravitational-wave driven CFS instability, which 
grows on a secular time-scale. Here we have found that while 
the actual mass-shedding-induced damping does not occur 
on dynamical time-scales (as in the Cowling approximation) 
it is still present on secular time-scales, so that a detailed 
comparison to the /-mode growth rate is required in order 
to determine the outcome of the /-mode instability. 

Taking into account that in the core collapse simula- 
tions by Dimmelmeier, Font & Miiller (2001, 2002a,b) the 
dominant quadrupole contribution had an effective gravita- 
tional wave amplitude that was roughly one to two orders of 
magnitude above the sensitivity curve of the advanced LIGO 
detector, it is likely that in a favourable detection event sev- 
eral pulsation modes of the proto-neutron star may be de- 
tected simultaneously. We have thus investigated the gravi- 
tational wave signals (in the Newtonian quadrupole approxi- 
mation) emitted by the linear modes F, Hi, ^f, and '^pi, for a 
slowly rotating neutron star model. In order to estimate the 
prospects for detectability of the signal by interferometric 
detectors, we have computed the (weakly detector depen- 
dent) characteristic signal frequency and amplitude for the 
theoretical power spectral density of VIRGO, LIGO I, and 
advanced LIGO. We infer that for a pulsation amplitude 
of the stellar central density of 5 per cent a gravitational 
wave signal from the investigated modes (with the excep- 
tion of the _ffi-mode) is measurable by current detectors, if 
the signal is integrated over a time of at least 20 ms and the 
source is located in our Galaxy. We have shown how scaling 
properties of the signal can be used to determine the neces- 
sary minimum stellar pulsation amplitude, maximum source 
distance, and minimum signal duration time required for a 
prospective detection. We have also discussed how effects of 
rotation can qualitatively influence empirical formulas used 
in gravitational wave asteroseismology and can thus help to 
constrain parameters of neutron star models in the case of 
a successful detection of gravitational waves from pulsation 
modes. A more complete study of the detectability of vari- 
ous modes for all equilibrium models of our sequences will 
appear elsewhere. 

In the relativistic core collapse simulations by Dim- 
melmeier, Font & Miiller (2001, 2002a,b) the quasi-periodic 
gravitational waves emitted around core bounce were found 
to have frequencies less than roughly 1.1 kHz. On the other 
hand, for nonrotating models constructed with a large sam- 
ple of different realistic EOSs, the frequencies of the fun- 
damental I = 2 quadrupole /-mode range from ~ 1.35 kHz 
for extremely stiff EOSs to ~ 3.6 kHz for extremely soft 
EOSs (see, e.g., Andersson & Kokkotas 1998). One of the 



main effects lowering the frequency of emitted gravitational 
waves is rotation. In our differentially rotating sequence of 
fixed rest mass, the nonrotating model has a fundamen- 
tal / — 2 /-mode frequency of ^ 1.6 kHz, which decreases 
to ~ 1.1 kHz at the largest rotation rate. Therefore, addi- 
tional effects must be taken into account in order to explain 
the wider range of frequencies found by Dimmelmeier, Font 
& Miiller (2001, 2002b). One such important factor is the 
high entropy in a proto-neutron star immediately after core 
bounce. The proto-neutron star is surrounded by a hot en- 
velope into which the pulsations of the interior can pene- 
trate, increasing the period of each pulsation (see Ferrari, 
Miniutti & Pons 2003a). Thus, a natural extension of our 
investigation will be to compute the frequencies of the var- 
ious axisymmetric modes for a realistic proto-neutron star 
structure. 
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